########################################################################
############ REPLICATION CODE FOR MCMURRY 2021 #########################
########################################################################

## Load packages
library(Matching)
library(plm)
library(lmtest)
library(spdep)
library(stargazer)
library(mvtnorm)
library(data.table)
library(AER)
library(ebal)
library(xtable)
library(sp)
library(rgdal)
library(dplyr)
library(multiwayvcov)
library(lfe)
library(fixest)
library(geosphere)
library(pastecs)

## Function to calculate Conley SEs
source("code/archived-code/deprecated-conley.R")

## Load data

# Administrative data

# Unique CADTs (for figure)
cadt_data <- read.csv("cadt_data_unique.csv", stringsAsFactors = F)
# Main barangay-level dataset
bgy <- read.csv("bgy_merged.csv", stringsAsFactors = F)

# Run matching code
source("bgy_matching.R")
# load("matched_samp_bgys.RData")

# Survey data

# Main survey dataset
all <- read.csv("community_survey.csv", stringsAsFactors = F)
# Conjoint results (Appendix only)
cjoint <- read.csv("survey_conjoint.csv", stringsAsFactors = F)
# Province and barangay maps (for survey sample maps in online appendix)
load("province_map.RData")
load("bgy_map.RData") 

# barangays with approved CADTs
table(bgy$Approved_CADT_NCIP_database[bgy$Approved_CADT_NCIP_database_munonly == 0])
# number of municipalities
length(unique(bgy$mun_2015cor[which(bgy$Approved_CADT_NCIP_database == 1)]))
# barangays with on-process CADTs or CADT applications
table(bgy$OP_or_FP[which(bgy$Approved_CADT_NCIP_database == 0 | bgy$Approved_CADT_NCIP_database_munonly == 1)])
# Identified as eligible - no other status
nrow(bgy[which(bgy$Identified_AD == 1 & bgy$Approved_CADT_NCIP_database == 0 & bgy$FP_direct == 0 & bgy$OP == 0), ])

#### TABLES ####

# Table 1: Land Titling (Continuous) and Indigenous Identification (Two-Way FE)

# Create panel analysis dataset with two periods
varyingvars_did <- c("ip2000", "ip2010",
                     "birth_reg2000", "birth_reg2010",
                     "leg_whipple2000", "leg_whipple2010",
                     "age2000", "age2010",
                     "hh_size2000", "hh_size2010",
                     "have_bhs2000", "have_bhs2010",
                     "have_highway2000", "have_highway2010",
                     "have_street_pattern2000", "have_street_pattern2010",
                     "same_mun2000", "same_mun2010",
                     "logpop2000", "logpop2010",
                     "cadt_prop2000", "cadt_prop2010",
                     "cadt_bin_update2000", "cadt_bin_update2010"
)
data_long <- reshape(bgy,
                     idvar = "GEOCODE2",
                     varying = varyingvars_did,
                     direction = "long",
                     sep = ""
                     )
# Drop barangays that are not present in both census years (changed borders or code)
data_long_nona <- data_long[!is.na(data_long$ip), ] 
data_long <- make.pbalanced(data_long_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))
# Create subsets for analysis
data_long_rural <- data_long[which(data_long$rural == 1 & data_long$reg_name != "CAR"), ]
data_long_matched <- data_long[data_long$GEOCODE2 %in% matched_samp_bgys, ]
data_long_universe <- data_long[data_long$universe_ncip_updated == 1 & data_long$reg_name != "CAR", ]
data_long_titled <- data_long[data_long$Approved_CADT_NCIP_database == 1 & data_long$Approved_CADT_NCIP_database_munonly == 0 & data_long$reg_name != "CAR", ]

# All rural
mod.ip.rural <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                   data = data_long_rural)
mod.ip.rural2 <- felm(ip ~ cadt_prop + same_mun + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                   data = data_long_rural)
# Matched
mod.ip.matched <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                     data = data_long_matched)
mod.ip.matched2 <- felm(ip ~ cadt_prop + same_mun + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                      data = data_long_matched)
# Universe
mod.ip.universe <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                       data = data_long_universe)
mod.ip.universe2 <- felm(ip ~ cadt_prop + same_mun + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                        data = data_long_universe)
# Titled
mod.ip.titled <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                        data = data_long_titled)
mod.ip.titled2 <- felm(ip ~ cadt_prop + same_mun + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                         data = data_long_titled)

# Create table
numobs <- c(nobs(mod.ip.rural) / 2,
            nobs(mod.ip.rural2) / 2,
            nobs(mod.ip.universe) / 2,
            nobs(mod.ip.universe2) / 2,
            nobs(mod.ip.matched) / 2,
            nobs(mod.ip.matched2) / 2,
            nobs(mod.ip.titled) / 2,
            nobs(mod.ip.titled2) / 2)
star <- stargazer(mod.ip.rural,
                  mod.ip.rural2,
                  mod.ip.universe,
                  mod.ip.universe2,
                  mod.ip.matched,
                  mod.ip.matched2,
                  mod.ip.titled,
                  mod.ip.titled2,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Indigenous Prop.",
                  font.size = "footnotesize",
                  column.labels = c("All Rural", "All Rural", "Eligible", "Eligible", "Matched", "Matched", "Titled", "Titled"),
                  covariate.labels = c("Titled Prop.",
                                       "Same Municipality",
                                       "Log. Population",
                                       "Mean Age",
                                       "Mean HH Size"),
                  add.lines = list(c("Observations", numobs)),
                  label = "tab:iddid",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling (Continuous) and Indigenous Identification")
cat(star, sep = '\n', 
    file = "id_continuous.tex"
    )

# coefficients in terms of standard deviations among untreated in each subset (for article text)

bgy$ip_delta <- bgy$ip_2010 - bgy$ip_2000

bgy_rural <- bgy[which(bgy$rural == 1 & bgy$reg_name != "CAR"), ]
bgy_universe <- bgy[bgy$universe_ncip_updated == 1 & bgy$reg_name != "CAR", ]
bgy_matched <- bgy[bgy$GEOCODE2 %in% matched_samp_bgys, ]
bgy_titled <- bgy[bgy$Approved_CADT_NCIP_database == 1 & bgy$Approved_CADT_NCIP_database_munonly == 0 & bgy$reg_name != "CAR", ]

summary(mod.ip.rural)$coefficients[1, 1] / sd(na.omit(bgy_rural$ip_delta[bgy_rural$cadt_bin_update2010 == 0]))
summary(mod.ip.matched)$coefficients[1, 1] / sd(na.omit(bgy_matched$ip_delta[bgy_matched$cadt_bin_update2010 == 0]))
summary(mod.ip.universe)$coefficients[1, 1] / sd(na.omit(bgy_universe$ip_delta[bgy_universe$cadt_bin_update2010 == 0]))
summary(mod.ip.titled)$coefficients[1, 1] / sd(na.omit(bgy_titled$ip_delta[bgy_titled$cadt_bin_update2010 == 0]))

# Table 2: Land Titling (Continuous) and Birth Registration (Two-Way FE)

# Create panel analysis dataset with four periods
varyingvars_panel <- c("birth_reg2000", "birth_reg2007", "birth_reg2010", "birth_reg2015",
                       "leg_whipple2000", "leg_whipple2007", "leg_whipple2010", "leg_whipple2015",
                       "have_bhs2000", "have_bhs2007", "have_bhs2010", "have_bhs2015",
                       "have_elem2000", "have_elem2007", "have_elem2010", "have_elem2015",
                       "have_hospital2000", "have_hospital2007", "have_hospital2010", "have_hospital2015",
                       "have_highway2000", "have_highway2007", "have_highway2010", "have_highway2015",
                       "have_hs2000", "have_hs2007", "have_hs2010", "have_hs2015",
                       "have_water_sys2000", "have_water_sys2007", "have_water_sys2010", "have_water_sys2015",
                       "have_street_pattern2000", "have_street_pattern2007", "have_street_pattern2010", "have_street_pattern2015",
                       "cadt_prop2000", "cadt_prop2007", "cadt_prop2010", "cadt_prop2015",
                       "cadt_bin_update2000", "cadt_bin_update2007", "cadt_bin_update2010", "cadt_bin_update2015",
                       "logpop2000", "logpop2007", "logpop2010", "logpop2015",
                       "hh_size2000", "hh_size2007", "hh_size2010", "hh_size2015",
                       "hs_grad2000", "hs_grad2007", "hs_grad2010", "hs_grad2015",
                       "age2000", "age2007", "age2010", "age2015")

data_long4 <- reshape(bgy,
                     idvar = "GEOCODE2",
                     varying = varyingvars_panel,
                     direction = "long",
                     sep = "")
# Drop barangays that are not present in both census years (changed borders or code)
data_long4_nona <- data_long4[!is.na(data_long4$birth_reg), ] 
data_long4 <- make.pbalanced(data_long4_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))

# Create subsets for analysis
data_long4_rural <- data_long4[which(data_long4$rural == 1 & data_long4$reg_name != "CAR"), ]
data_long4_matched <- data_long4[data_long4$GEOCODE2 %in% matched_samp_bgys, ]
data_long4_universe <- data_long4[data_long4$universe_ncip_updated == 1 & data_long4$reg_name != "CAR", ]
data_long4_titled <- data_long4[data_long4$Approved_CADT_NCIP_database == 1 & data_long4$Approved_CADT_NCIP_database_munonly == 0 & data_long4$reg_name != "CAR", ]

# barangays missing barangay-level PG variables
data_long4_nona_controls <- data_long4[!is.na(data_long4$have_bhs) & !is.na(data_long4$have_street_pattern) & !is.na(data_long4$have_highway), ]
data_long4_controls <- make.pbalanced(data_long4_nona_controls, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))

# Create subsets for analysis
data_long4_controls_rural <- data_long4_controls[data_long4_controls$rural == 1 & data_long4_controls$reg_name != "CAR", ]
data_long4_controls_matched <- data_long4_controls[data_long4_controls$GEOCODE2 %in% matched_samp_bgys, ]
data_long4_controls_universe <- data_long4_controls[data_long4_controls$universe_ncip_updated == 1 & data_long4_controls$reg_name != "CAR", ]
data_long4_controls_titled <- data_long4_controls[data_long4_controls$Approved_CADT_NCIP_database == 1 & data_long4_controls$Approved_CADT_NCIP_database_munonly == 0 & data_long4_controls$reg_name != "CAR", ]

# All rural
mod.breg.rural <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                       data = data_long4_rural)
mod.breg.rural2 <- felm(birth_reg ~ cadt_prop + have_bhs + have_street_pattern + have_highway + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                        data = data_long4_controls_rural)
# Matched
mod.breg.matched <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                         data = data_long4_matched)
mod.breg.matched2 <- felm(birth_reg ~ cadt_prop + have_bhs + have_street_pattern + have_highway + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                          data = data_long4_controls_matched)
# Universe
mod.breg.universe <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                          data = data_long4_universe)
mod.breg.universe2 <- felm(birth_reg ~ cadt_prop + have_bhs + have_street_pattern + have_highway + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                           data = data_long4_controls_universe)
# Titled
mod.breg.titled <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                        data = data_long4_titled)
mod.breg.titled2 <- felm(birth_reg ~ cadt_prop + have_bhs + have_street_pattern + have_highway + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                         data = data_long4_controls_titled)

# Create table
numobs <- c(nobs(mod.breg.rural) / 4,
            nobs(mod.breg.rural2) / 4,
            nobs(mod.breg.universe) / 4,
            nobs(mod.breg.universe2) / 4,
            nobs(mod.breg.matched) / 4,
            nobs(mod.breg.matched2) / 4,
            nobs(mod.breg.titled) / 4,
            nobs(mod.breg.titled2) / 4)
star <- stargazer(mod.breg.rural,
                  mod.breg.rural2,
                  mod.breg.universe,
                  mod.breg.universe2,
                  mod.breg.matched,
                  mod.breg.matched2,
                  mod.breg.titled,
                  mod.breg.titled2,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Birth Registration",
                  font.size = "footnotesize",
                  column.labels = c("All Rural", "All Rural", "Eligible", "Eligible", "Matched", "Matched", "Titled", "Titled"),
                  covariate.labels = c("Titled Prop.",
                                       "Bgy. Health Center",
                                       "Street Pattern",
                                       "Highway Access",
                                       "Log Population",
                                       "Mean Age",
                                       "Mean HH Size"),
                  add.lines = list(c("Observations", numobs)),
                  label = "tab:bregpanel",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling (Continuous) and Birth Registration (Two-Way FE)")
cat(star, sep = '\n', 
    file = "breg_continuous.tex"
    )

# coefficients in terms of standard deviations among untreated in each subset (for article text)

bgy$birth_reg_delta <- bgy$birth_reg_2010 - bgy$birth_reg_2000

bgy_rural <- bgy[which(bgy$rural == 1 & bgy$reg_name != "CAR"), ]
bgy_universe <- bgy[which(bgy$universe_ncip_updated == 1 & bgy$reg_name != "CAR"), ]
bgy_matched <- bgy[which(bgy$GEOCODE2 %in% matched_samp_bgys), ]
bgy_titled <- bgy[which(bgy$Approved_CADT_NCIP_database == 1 & bgy$Approved_CADT_NCIP_database_munonly == 0 & bgy$reg_name != "CAR"), ]

# coefficients in terms of standard deviations among untreated

summary(mod.breg.rural)$coefficients[1, 1] / sd(na.omit(bgy_rural$birth_reg_delta[bgy_rural$cadt_bin_update2015 == 0]))
summary(mod.breg.matched)$coefficients[1, 1] / sd(na.omit(bgy_matched$birth_reg_delta[bgy_matched$cadt_bin_update2015 == 0]))
summary(mod.breg.universe)$coefficients[1, 1] / sd(na.omit(bgy_universe$birth_reg_delta[bgy_universe$cadt_bin_update2015 == 0]))
summary(mod.breg.titled)$coefficients[1, 1] / sd(na.omit(bgy_titled$birth_reg_delta[bgy_titled$cadt_bin_update2015 == 0]))

# in terms of baseline differences in birth registration

# non-indigenous majority vs majority indigenous rural barangays
maj_diff_rural <- mean(bgy_rural$birth_reg_2000[which(bgy_rural$ip_2000 < 0.5)], na.rm = T) - mean(bgy_rural$birth_reg_2000[which(bgy_rural$ip_2000 > 0.5)], na.rm = T)

summary(mod.breg.rural)$coefficients[1, 1] / maj_diff_rural
summary(mod.breg.matched)$coefficients[1, 1] / maj_diff_rural
summary(mod.breg.universe)$coefficients[1, 1] / maj_diff_rural
summary(mod.breg.titled)$coefficients[1, 1] / maj_diff_rural

# Pre-treatment trends among barangays titled after 2010 (for article text)
bgy$cadt_approval_year_first_update <- bgy$cadt_approval_year_first
bgy$cadt_approval_year_first_update[bgy$Approved_CADT_NCIP_database_munonly == 1] <- NA
bgy$ip_delta <- bgy$ip2010 - bgy$ip2000
# subset to those that received approval post-2010
bgy_after <- bgy[bgy$cadt_approval_year_first_update > 2010 & !is.na(bgy$cadt_approval_year_first_update) & bgy$reg_name != "CAR", ]
bgy_after$cadt_approval_year_first_update[bgy_after$cadt_approval_year_first_update == "Inf"] <- NA
# comparing 2000-2010 change for barangays that got CADTs 2010 to 2014
# vs those who got 2014 - 2018
t.test(bgy_after$ip_delta[bgy_after$cadt_approval_year_first_update <= 2014], bgy_after$ip_delta[bgy_after$cadt_approval_year_first_update > 2014])

# Table 3: Prime Treatment and Top Ranking of Identity Attributes

# subset to only identity prime and control (see explanation in Appendix)
all_nomat <- all[all$ipra_mat == 0, ]

# demeaned covariates
all_nomat$edu_hs_dm <- all_nomat$edu_hs - mean(all_nomat$edu_hs, na.rm = T)
all_nomat$rel_evangelical_dm <- all_nomat$rel_evangelical - mean(all_nomat$rel_evangelical, na.rm = T)

all_nomat$tribe_top_placement_dm <- all_nomat$tribe_top_placement - mean(all_nomat$tribe_top_placement, na.rm = T)
all_nomat$nationality_top_placement_dm <- all_nomat$nationality_top_placement - mean(all_nomat$nationality_top_placement, na.rm = T)
all_nomat$religion_top_placement_dm <- all_nomat$religion_top_placement - mean(all_nomat$religion_top_placement, na.rm = T)
all_nomat$gender_top_placement_dm <- all_nomat$gender_top_placement - mean(all_nomat$gender_top_placement, na.rm = T)

all_nomat$tribe_placement_dm <- all_nomat$tribe_placement - mean(all_nomat$tribe_placement, na.rm = T)
all_nomat$nationality_placement_dm <- all_nomat$nationality_placement - mean(all_nomat$nationality_placement, na.rm = T)
all_nomat$religion_placement_dm <- all_nomat$religion_placement - mean(all_nomat$religion_placement, na.rm = T)
all_nomat$gender_placement_dm <- all_nomat$gender_placement - mean(all_nomat$gender_placement, na.rm = T)

# without covariates
mod.tribetop.id <- lm(tribe_top ~ ipra_any, data = all_nomat)
mod.nattop.id <- lm(nationality_top ~ ipra_any, data = all_nomat)
mod.gendertop.id <- lm(gender_top ~ ipra_any, data = all_nomat)
mod.religiontop.id <- lm(religion_top ~ ipra_any, data = all_nomat)

# covariate adjusted
mod.tribetop.id.cov <- lm(tribe_top ~ ipra_any*(edu_hs_dm + rel_evangelical_dm + tribe_placement_dm), data = all_nomat)
mod.nattop.id.cov <- lm(nationality_top ~ ipra_any*(edu_hs_dm + rel_evangelical_dm + nationality_placement_dm), data = all_nomat)
mod.gendertop.id.cov <- lm(gender_top ~ ipra_any*(edu_hs_dm + rel_evangelical_dm + gender_placement_dm), data = all_nomat)
mod.religiontop.id.cov <- lm(religion_top ~ ipra_any*(edu_hs_dm + rel_evangelical_dm + religion_placement_dm), data = all_nomat)

# Create table
stargazer(mod.tribetop.id,
          mod.tribetop.id.cov,
          mod.nattop.id,
          mod.nattop.id.cov,
          mod.gendertop.id,
          mod.gendertop.id.cov,
          mod.religiontop.id,
          mod.religiontop.id.cov,
          font.size = "footnotesize",
          title = "Prime Treatment and Top Ranking of Identity Attributes",
          covariate.labels = c("IPRA Prime"),
          dep.var.labels = c("Tribe Top", "Nationality Top", "Gender Top", "Religion Top"),
          omit = c(":", "edu", "rel", "placement", "Constant"),
          omit.stat = c("f", "ser"),
          add.lines = list(c("Covariate Adjustment", rep(c("N", "Y"), 8))),
          type = "latex",
          label = "tab:alltopidonly",
          out = "alltopidonly.tex"
)

# Effect of prime on belief government conducted survey (for article text)
mod.conducted.gov.id <- lm(conducted_gov ~ ipra_any, data = all_nomat)
summary(mod.conducted.gov.id)

#### FIGURES ####

# Figure 1: Cumulative number of Certificates of Ancestral Domain Title (CADT) approved per year

pdf("cadt_by_year.pdf")
barplot(cumsum(table(cadt_data$approved_year)), xlab = "Year", ylab = "Total approved CADTs", main = "Total Approved CADTs by Year")
dev.off()

# Figure 2: Land area covered by approved Certificates of Ancestral Domain Title (CADT) as of October 2016 and eligible areas

## SHARING OF SHAPEFILES USED TO REPRODUCE THIS FIGURE REQUIRES PERMISSION FROM THE PHILIPPINES'
## NATIONAL COMMISSION ON INDIGENOUS PEOPLES, PLEASE CONTACT AUTHOR FOR MORE INFORMATION.

#### APPENDIX TABLES ####

# Table A1: Pre-Treatment Covariate Balance Among Eligible Barangays (2018)

bgy_universe <- bgy[which(bgy$universe_ncip_updated == 1 & bgy$reg_name != "CAR"), ]

covars <- c("ip_delta", "GEOCODE2", "Prov_Code_num", "Mun_Code_num", 
            "ip_2000", "birth_reg_2000", "ethfrac_2000", "area_sqkm", "elev_mean", "elev_sd",
            "slope_mean", "soil_index", "mindep_bin", "log_ncip_dist", "log_coast_dist", "log_road_dist", "religion_catholic_2000",
            "have_street_pattern_2000", "have_highway_2000", "have_church_2000",
            "have_market_2000", "have_elem_2000", "have_bhs_2000", "have_water_sys_2000",
            "leg_pretrend", "hs_grad_pretrend", "elem_grad_pretrend", "logpop_pretrend", "strong_mat_pretrend")

exclude <- c("ip_delta", "GEOCODE2", "Prov_Code_num", "Mun_Code_num")
baltab2018_universe <- baltest.collect(MatchBalance(formula(paste("cadt_bin_update2018 ~", paste0(covars[-which(covars %in% exclude)], collapse = "+"), sep = " " )), data = bgy_universe), var.names = covars[-which(covars %in% exclude)], after = FALSE)
baltab2018_universe_sub <- baltab2018_universe[, c(1, 2, 6)]
colnames(baltab2018_universe_sub) <- c("Mean Titled", "Mean Untitled", "T p-val")
rownames(baltab2018_universe_sub) <- c("Indigenous Prop. 2000",
                          "Birth Registration 2000",
                          "Ethnic Frac. 2000",
                          "Area (sq. km)",
                          "Elevation Mean",
                          "Elevation St. Dev",
                          "Slope Mean",
                          "Soil Quality Index",
                          "Mineral Deposits",
                          "Log. NCIP Dist.",
                          "Log. Coast Dist.",
                          "Log. Road Dist.",
                          "Catholic Prop. 2000",
                          "Street Pattern 2000",
                          "Highway Access 2000",
                          "Church 2000",
                          "Market 2000",
                          "Elementary 2000",
                          "Bgy. Health Ctr 2000",
                          "Water System 2000",
                          "Legibility $\\Delta$ 1990-2000",
                          "HS Grad. $\\Delta$ 1990-2000",
                          "Elem Grad. $\\Delta$ 1990-2000",
                          "Log Pop $\\Delta$ 1990-2000",
                          "Housing Quality $\\Delta$ 1990-2000"
)
baltab2018_universex <- xtable(baltab2018_universe_sub, 
                               caption = "Pre-Treatment Covariate Balance Among Eligible Barangays (Titled in 2018 vs. Untitled in 2018)",
                               digits = 3,
                               label = "tab:unibalance")
print(baltab2018_universex,
      sanitize.text.function = function(x){x},
      type = "latex",
      file = "universe_balance.tex",
      caption.placement = "top"
)

# Table A.2: Pre-treatment Covariate Balance Among Barangays Titled pre-2010 and Comparison Groups

bgy_universe <- bgy[which(bgy$universe_ncip_updated == 1 & bgy$reg_name != "CAR"), ]
bgy_titled <- bgy[which(bgy$Approved_CADT_NCIP_database == 1 & bgy$Approved_CADT_NCIP_database_munonly == 0 & bgy$reg_name != "CAR"), ]
bgy_matched <- bgy[which(bgy$GEOCODE2 %in% matched_samp_bgys), ]

covars <- c("ip_delta", "GEOCODE2", "Prov_Code_num", "Mun_Code_num", 
            "ip_2000", "birth_reg_2000", "ethfrac_2000", "area_sqkm", "elev_mean", "elev_sd",
            "slope_mean", "soil_index", "mindep_bin", "log_ncip_dist", "log_coast_dist", "log_road_dist", "religion_catholic_2000",
            "have_street_pattern_2000", "have_highway_2000", "have_church_2000",
            "have_market_2000", "have_elem_2000", "have_bhs_2000", "have_water_sys_2000",
            "leg_pretrend", "hs_grad_pretrend", "elem_grad_pretrend", "logpop_pretrend", "strong_mat_pretrend")
exclude <- c("ip_delta", "GEOCODE2", "Prov_Code_num", "Mun_Code_num")

baltab2010_universe <- baltest.collect(MatchBalance(formula(paste("cadt_bin_update2010 ~", paste0(covars[-which(covars %in% exclude)], collapse = "+"), sep = " " )), data = bgy_universe, print.level = 0), var.names = covars[-which(covars %in% exclude)], after = FALSE)
baltab2010_universe_sub <- baltab2010_universe[, c(1, 2, 6)]

baltab2010_titled <- baltest.collect(MatchBalance(formula(paste("cadt_bin_update2010 ~", paste0(covars[-which(covars %in% exclude)], collapse = "+"), sep = " " )), data = bgy_titled, print.level = 0), var.names = covars[-which(covars %in% exclude)], after = FALSE)
baltab2010_titled_sub <- baltab2010_titled[, c(1, 2, 6)]

baltab2010_matched <- baltest.collect(MatchBalance(formula(paste("cadt_bin_update2010 ~", paste0(covars[-which(covars %in% exclude)], collapse = "+"), sep = " " )), data = bgy_matched, print.level = 0), var.names = covars[-which(covars %in% exclude)], after = FALSE)
baltab2010_matched_sub <- baltab2010_matched[, c(1, 2, 6)]

baltab_all <- cbind(baltab2010_universe[, c(1, 2)],
                    baltab2010_titled[, 2],
                    baltab2010_matched[, 2],
                    baltab2010_universe[, 6],
                    baltab2010_titled[, 6],
                    baltab2010_matched[, 6]
)
colnames(baltab_all) <- c("Titled pre 2010",
                          "Eligible Untitled pre 2010",
                          "Titled post 2010",
                          "Matched Controls",
                          "T p-val: Titled pre vs Eligible",
                          "T p-val: Titled post vs Titled pre",
                          "T p-val: Titled pre vs. Matched")
rownames(baltab_all) <- c("Indigenous Prop. 2000",
                          "Birth Registration 2000",
                          "Ethnic Frac. 2000",
                          "Area (sq. km)",
                          "Elevation Mean",
                          "Elevation St. Dev",
                          "Slope Mean",
                          "Soil Quality Index",
                          "Mineral Deposits",
                          "Log. NCIP Dist.",
                          "Log. Coast Dist.",
                          "Log. Road Dist.",
                          "Catholic Prop. 2000",
                          "Street Pattern 2000",
                          "Highway Access 2000",
                          "Church 2000",
                          "Market 2000",
                          "Elementary 2000",
                          "Bgy. Health Ctr 2000",
                          "Water System 2000",
                          "Legibility $\\Delta$ 1990-2000",
                          "HS Grad. $\\Delta$ 1990-2000",
                          "Elem Grad. $\\Delta$ 1990-2000",
                          "Log Pop $\\Delta$ 1990-2000",
                          "Housing Quality $\\Delta$ 1990-2000"
)
baltab_allx <- xtable(baltab_all, 
                      caption = "Pre-Treatment Covariate Balance Among Barangays Titled pre-2010 and Comparison Groups",
                      digits = 3,
                      align = c("p{5cm}",
                                rep("p{1.75cm}", 3),
                                "p{1.75cm}|",
                                rep("p{1.75cm}", 3)),
                      label = "tab:allbalance")
print(baltab_allx,
      sanitize.text.function = function(x){x},
      type = "latex",
      file = "all_balance.tex",
      caption.placement = "top"
)

# Table A.3: Summary Statistics (Study Barangays)

vars <- c("cadt_prop2010", "cadt_bin_update2010", "cadt_bin_update2018", "universe_ncip_updated",
          "ip_2000", "ip_2010", "birth_reg_2000", "birth_reg_2007", "birth_reg_2010", "birth_reg_2015")
bgy_rural <- bgy[which(bgy$rural == 1 & bgy$reg_name != "CAR"), ]
bgy_rural_sub <- bgy_rural[, vars]

deletecols <- c("nbr.null", "range", "SE.mean", "CI.mean.0.95", "var", "coef.var", "sum")
sumtab <- t(stat.desc(bgy_rural_sub))
sumtabsub <- sumtab[, -which(colnames(sumtab) %in% deletecols)]
mdat <- matrix(c(rep(0, nrow(sumtabsub)),
                 rep(0, nrow(sumtabsub)),
                 rep(0, nrow(sumtabsub)),
                 rep(0, nrow(sumtabsub)),
                 rep(0, nrow(sumtabsub)),
                 rep(2, nrow(sumtabsub)),
                 rep(2, nrow(sumtabsub)),
                 rep(2, nrow(sumtabsub))),
               nrow = nrow(sumtabsub),
               ncol = ncol(sumtabsub) + 1,
               byrow = F)

sumtabx <- xtable(sumtabsub, 
                  caption = "Summary Statistics (Study Barangays)",
                  label = "tab:sumstats",
                  digits = mdat)

rownames(sumtabx) <- c("Titled Proportion 2010",
                       "Titled Binary 2010",
                       "Titled Binary 2018",
                       "Eligible Universe",
                       "Indigenous Pop. 2000",
                       "Indigenous Pop. 2010",
                       "Birth Registration 2000",
                       "Birth Registration 2007",
                       "Birth Registration 2010",
                       "Birth Registration 2015")
colnames(sumtabx) <- c("Num. Vals", "Num. NA", "Min", "Max", "Median", "Mean", "Std. Dev")

print(sumtabx,
      type = "latex",
      file = "summary_stats.tex",
      caption.placement = "top",
      font.size = "footnotesize")

# Table A.4: Land Titling and Indigenous Identity (Quasi-RDD)

bgy$cadt_approval_year_first_update <- bgy$cadt_approval_year_first
bgy$cadt_approval_year_first_update[bgy$Approved_CADT_NCIP_database_munonly == 1] <- NA

# Change from 2000 to 2010 in DVs
bgy$ip_delta <- bgy$ip_2010 - bgy$ip_2000

# Create subsets
bgy_rural <- bgy[which(bgy$rural == 1 & bgy$reg_name != "CAR"), ]
bgy_universe <- bgy[which(bgy$universe_ncip_updated == 1 & bgy$reg_name != "CAR"), ]
bgy_matched <- bgy[which(bgy$GEOCODE2 %in% matched_samp_bgys), ]
bgy_titled <- bgy[which(bgy$Approved_CADT_NCIP_database == 1 & bgy$Approved_CADT_NCIP_database_munonly == 0 & bgy$reg_name != "CAR"), ]

# Distance from titling year to 2010 census
bgy_titled$census_distance <- bgy_titled$cadt_approval_year_first_update - 2010

# All barangays (8-year bandwidth)
mod.quasi.rdd.all <- lm(ip_delta ~ cadt_bin_update2010*census_distance, data = bgy_titled[bgy_titled$census_distance %in% c(-8:8), ])

# 4-year bandwidth
mod.quasi.rdd.all4 <- lm(ip_delta ~ cadt_bin_update2010*census_distance, data = bgy_titled[bgy_titled$census_distance %in% c(-4:4), ])

stargazer(mod.quasi.rdd.all4,
          mod.quasi.rdd.all,
          covariate.labels = c("Titled",
                               "Census Distance",
                               "Census Distance x Titled",
                               "Constant"),
          omit.stat = c("ser", "f"),
          dep.var.labels = "Indigenous Prop. $\\Delta$ 2000-2010",
          add.lines = list(c("Bandwidth", 4, 8)),
          title = "Land Titling and Indigenous Identity (Quasi-RDD)",
          label = "tab:idrddbin",
          out = "id_rdd_bin.tex"
          )

# Table A.5: Land Titling and Indigenous Identification (First-difference with filing year)

# Restrict to barangays associated with a CADT not missing a filing year
bgy_filed <- bgy[!is.na(bgy$filed_year) & bgy$reg_name != "CAR", ]

# control for filed year as a continuous variable
ip.delta.filed <- lm(ip_delta ~ cadt_prop2010 + filed_year, data = bgy_filed)

# filing year fixed effects
ip.delta.filed2 <- lm(ip_delta ~ cadt_prop2010 + as.factor(filed_year), data = bgy_filed)

# repeat with binary titling variable
ipbin.delta.filed <- lm(ip_delta ~ cadt_bin_update2010 + filed_year, data = bgy_filed)

ipbin.delta.filed2 <- lm(ip_delta ~ cadt_bin_update2010 + as.factor(filed_year), data = bgy_filed)

# Create table
star <- stargazer(ip.delta.filed,
                  ip.delta.filed2,
                  ipbin.delta.filed,
                  ipbin.delta.filed2,
                  omit = c("as.factor",
                           "Constant"),
                  omit.stat = c("ser", "f"),
                  dep.var.labels = "Indigenous Prop. $\\Delta$ 2000-2010",
                  covariate.labels = c("Titled Prop.",
                                       "Titled (Binary)",
                                       "Filed Year"),
                  label = "tab:idfile",
                  add.lines = list(c("Filed Year FE", "N", "Y", "N", "Y")),
                  title = "Land Titling and Indigenous Identification (First-difference with filing year)")
cat(star, 
    sep = "\n", 
    file = "id_filing.tex"
    )

# Table A.6: Land Titling (Binary) and Indigenous Identification

# All rural
mod.ipbin.rural <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                        data = data_long_rural)
mod.ipbin.rural2 <- felm(ip ~ cadt_bin_update + same_mun + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                         data = data_long_rural)
# Matched
mod.ipbin.matched <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                          data = data_long_matched)
mod.ipbin.matched2 <- felm(ip ~ cadt_bin_update + same_mun + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                           data = data_long_matched)
# Universe
mod.ipbin.universe <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                           data = data_long_universe)
mod.ipbin.universe2 <- felm(ip ~ cadt_bin_update + same_mun + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                            data = data_long_universe)
# Titled
mod.ipbin.titled <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                         data = data_long_titled)
mod.ipbin.titled2 <- felm(ip ~ cadt_bin_update + same_mun + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                          data = data_long_titled)

# Create table
numobs <- c(nobs(mod.ipbin.rural) / 2,
            nobs(mod.ipbin.rural2) / 2,
            nobs(mod.ipbin.universe) / 2,
            nobs(mod.ipbin.universe2) / 2,
            nobs(mod.ipbin.matched) / 2,
            nobs(mod.ipbin.matched2) / 2,
            nobs(mod.ipbin.titled) / 2,
            nobs(mod.ipbin.titled2) / 2)
star <- stargazer(mod.ipbin.rural,
                  mod.ipbin.rural2,
                  mod.ipbin.universe,
                  mod.ipbin.universe2,
                  mod.ipbin.matched,
                  mod.ipbin.matched2,
                  mod.ipbin.titled,
                  mod.ipbin.titled2,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Indigenous Prop.",
                  font.size = "footnotesize",
                  column.labels = rep(c("All Rural", "Eligible", "Matched", "Titled"), each = 2),
                  covariate.labels = c("Titled (Binary)",
                                       "Same Municipality",
                                       "Log. Population",
                                       "Mean Age",
                                       "Mean HH Size"),
                  add.lines = list(c("Observations", numobs)),
                  label = "tab:iddidbin",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling (Binary) and Indigenous Identification")
cat(star, sep = '\n', 
    file = "id_bin.tex"
    )

# Table A.7: Land Titling and Population Change 2000-2010

# All rural
mod.pop.rural <- felm(logpop ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                      data = data_long_rural)
# Matched
mod.pop.matched <- felm(logpop ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                        data = data_long_matched)
# Universe
mod.pop.universe <- felm(logpop ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                         data = data_long_universe)
# Titled
mod.pop.titled <- felm(logpop ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                       data = data_long_titled)

# Create table
numobs <- c(nobs(mod.pop.rural) / 2,
            nobs(mod.pop.universe) / 2,
            nobs(mod.pop.matched) / 2,
            nobs(mod.pop.titled) / 2)

star <- stargazer(mod.pop.rural,
                  mod.pop.universe,
                  mod.pop.matched,
                  mod.pop.titled,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Logged Population",
                  font.size = "footnotesize",
                  column.labels = c("All Rural", "Eligible", "Matched", "Titled"),
                  covariate.labels = "Titled Prop.",
                  add.lines = list(c("Observations", numobs)),
                  label = "tab:popchangedid",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling and Population")
cat(star, sep = '\n',
    file = "popchange.tex"
    )

# Table A.8: Land Titling and Recent Migration 2000-2010

# All rural
mod.mig.rural <- felm(same_mun ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                      data = data_long_rural)
# Matched
mod.mig.matched <- felm(same_mun ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                        data = data_long_matched)
# Universe
mod.mig.universe <- felm(same_mun ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                         data = data_long_universe)
# Titled
mod.mig.titled <- felm(same_mun ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                       data = data_long_titled)

# Create table
numobs <- c(nobs(mod.mig.rural) / 2,
            nobs(mod.mig.universe) / 2,
            nobs(mod.mig.matched) / 2,
            nobs(mod.mig.titled) / 2)

star <- stargazer(mod.mig.rural,
                  mod.mig.universe,
                  mod.mig.matched,
                  mod.mig.titled,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Same Municipality",
                  font.size = "footnotesize",
                  column.labels = c("All Rural", "Eligible", "Matched", "Titled"),
                  covariate.labels = "Titled Prop.",
                  add.lines = list(c("Observations", numobs)),
                  label = "tab:samemundid",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling and Recent Migration")
cat(star, sep = '\n',
   file = "samemun.tex"
)

# Table A.9: Land Titling and Indigenous Identification (First Differences with Province and Region Fixed Effects)

bgy$ip_delta <- bgy$ip2010 - bgy$ip_2000

bgy_rural <- bgy[which(bgy$rural == 1 & bgy$reg_name != "CAR"), ]
bgy_universe <- bgy[which(bgy$universe_ncip_updated == 1 & bgy$reg_name != "CAR"), ]
bgy_matched <- bgy[which(bgy$GEOCODE2 %in% matched_samp_bgys), ]
bgy_titled <- bgy[which(bgy$Approved_CADT_NCIP_database == 1 & bgy$Approved_CADT_NCIP_database_munonly == 0 & bgy$reg_name != "CAR"), ]

# Province FE
ip.delta.prov <- lm(ip_delta ~ cadt_prop2010 + I(prov_code_2000), data = bgy_rural)

ip.delta.prov.matched <- lm(ip_delta ~ cadt_prop2010 + I(prov_code_2000), data = bgy_matched)

ip.delta.prov.universe <- lm(ip_delta ~ cadt_prop2010 + I(prov_code_2000), data = bgy_universe)

ip.delta.prov.titled <- lm(ip_delta ~ cadt_prop2010 + I(prov_code_2000), data = bgy_titled)

# Region FE
ip.delta.reg <- lm(ip_delta ~ cadt_prop2010 + I(reg_name), data = bgy_rural)

ip.delta.reg.matched <- lm(ip_delta ~ cadt_prop2010 + I(reg_name), data = bgy_matched)

ip.delta.reg.universe <- lm(ip_delta ~ cadt_prop2010 + I(reg_name), data = bgy_universe)

ip.delta.reg.titled <- lm(ip_delta ~ cadt_prop2010 + I(reg_name), data = bgy_titled)

# Create table
star <- stargazer(ip.delta.prov,
                  ip.delta.prov.universe,
                  ip.delta.prov.matched,
                  ip.delta.prov.titled,
                  ip.delta.reg,
                  ip.delta.reg.universe,
                  ip.delta.reg.matched,
                  ip.delta.reg.titled,
                  omit = c("I\\(", "Constant"),
                  omit.stat = c("rsq", "f", "ser"),
                  type = "latex",
                  font.size = "footnotesize",
                  label = "tab:iddidprovreg",
                  column.labels = rep(c("All Rural", "Eligible", "Matched", "Titled"), 2),
                  dep.var.labels = "Indigenous Prop. $\\Delta$ 2000-2010",
                  covariate.labels = "Titled Prop.",
                  add.lines = list(c("Prov FE", c(rep("Y", 4), rep("N", 4))),
                                   c("Reg FE", c(rep("N", 4), rep("Y", 4)))
                                   ),
                  title = "Land Titling and Indigenous Identification (First-difference with Province and Region Fixed Effects)")
cat(star, sep = '\n', 
    file = "id_provreg.tex"
    )

# Table A.10: Land Titling and Indigenous Identification (Conley SE)

mod.ip.rural.conley <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | lat + lon, data = data_long_rural, keepCX = TRUE)
SE <- ConleySEs(reg = mod.ip.rural.conley,
                unit = "GEOCODE2",
                time = "time",
                lat = "lat", lon = "lon",
                dist_fn = "SH", dist_cutoff = 500,
                lag_cutoff = 5,
                cores = 1,
                verbose = FALSE)
ses.ip.rural <- sapply(SE, function(x) diag(sqrt(x))) %>% round(3)

mod.ip.matched.conley <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | lat + lon, data = data_long_matched, keepCX = TRUE)
SE_matched <- ConleySEs(reg = mod.ip.matched.conley,
                        unit = "GEOCODE2",
                        time = "time",
                        lat = "lat", lon = "lon",
                        dist_fn = "SH", dist_cutoff = 500,
                        lag_cutoff = 5,
                        cores = 1,
                        verbose = FALSE)
ses.ip.matched <- sapply(SE_matched, function(x) diag(sqrt(x))) %>% round(3)

mod.ip.universe.conley <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | lat + lon, data = data_long_universe, keepCX = TRUE)
SE_universe <- ConleySEs(reg = mod.ip.universe.conley,
                         unit = "GEOCODE2",
                         time = "time",
                         lat = "lat", lon = "lon",
                         dist_fn = "SH", dist_cutoff = 500,
                         lag_cutoff = 5,
                         cores = 1,
                         verbose = FALSE)
ses.ip.universe <- sapply(SE_universe, function(x) diag(sqrt(x))) %>% round(3)

mod.ip.titled.conley <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | lat + lon, data = data_long_titled, keepCX = TRUE)
SE_titled <- ConleySEs(reg = mod.ip.titled.conley,
                       unit = "GEOCODE2",
                       time = "time",
                       lat = "lat", lon = "lon",
                       dist_fn = "SH", dist_cutoff = 500,
                       lag_cutoff = 5,
                       cores = 1,
                       verbose = FALSE)
ses.ip.titled <- sapply(SE_titled, function(x) diag(sqrt(x))) %>% round(3)

# Create table
ses.all <- c(ses.ip.rural[3], ses.ip.matched[3], ses.ip.universe[3], ses.ip.titled[3])

numobs <- c(nobs(mod.ip.rural.conley) / 2,
            nobs(mod.ip.universe.conley) / 2,
            nobs(mod.ip.matched.conley) / 2,
            nobs(mod.ip.titled.conley) / 2)

rsqs <- round(c(summary(mod.ip.rural.conley)$r.squared,
                summary(mod.ip.universe.conley)$r.squared,
                summary(mod.ip.matched.conley)$r.squared,
                summary(mod.ip.titled.conley)$r.squared
), digits = 3)

star <- stargazer(mod.ip.rural.conley,
                  mod.ip.universe.conley,
                  mod.ip.matched.conley,
                  mod.ip.titled.conley,
                  se = ses.all,
                  omit.stat = c("rsq", "f", "ser", "n", "adj.rsq"),
                  type = "latex",
                  font.size = "footnotesize",
                  label = "tab:iddidconley",
                  column.labels = c("All Rural", "Eligible", "Matched", "Titled"),
                  dep.var.labels = "Indigenous Prop.",
                  covariate.labels = "Titled Prop.",
                  add.lines = list(c("Observations", numobs),
                                   c("$R^2$", rsqs)),
                  title = "Land Titling and Indigenous Identification (Conley SE)")
cat(star, sep = '\n',
   file = "id_did_conley.tex"
    )

# Table A.11: Land Titling and Indigenous Identification - Source Agreement

bgy_universe <- bgy[bgy$universe_ncip_updated == 1 & bgy$reg_name != "CAR", ]

# map says over 50% overlap and list says 1
confirmed_maj_2016 <- bgy_universe$GEOCODE2[bgy_universe$cadt_prop2016 > 0.5 & bgy_universe$cadt_bin_update2016 == 1]
# map says total overlap and list says 1 (in 2016)
confirmed_total_2016 <- bgy_universe$GEOCODE2[bgy_universe$cadt_prop2016 == 1 & bgy_universe$cadt_bin_update2016 == 1]
# map says zero overlap and list says zero overlap (in 2016)
confirmed_none_2016 <- bgy_universe$GEOCODE2[bgy_universe$cadt_prop2016 == 0 & bgy_universe$cadt_bin_update2016 == 0]

bgy_universe_none_total_2016 <- c(confirmed_total_2016, confirmed_none_2016)
bgy_universe_none_maj_2016 <- c(confirmed_maj_2016, confirmed_none_2016)

# Universe
# majority overlap - continuous
mod.ip.universe.confirmed.maj <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                                      data = data_long_universe[data_long_universe$GEOCODE2 %in% bgy_universe_none_maj_2016, ])
# majority overlap - binary
mod.ipbin.universe.confirmed.maj <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                         data = data_long_universe[data_long_universe$GEOCODE2 %in% bgy_universe_none_maj_2016, ])

# total overlap (binary only)
mod.ipbin.universe.confirmed.total <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                           data = data_long_universe[data_long_universe$GEOCODE2 %in% bgy_universe_none_total_2016, ])

# Titled
# majority overlap - continuous
mod.ip.titled.confirmed.maj <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                                    data = data_long_titled[data_long_titled$GEOCODE2 %in% bgy_universe_none_maj_2016, ])

# majority overlap - binary
mod.ipbin.titled.confirmed.maj <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                       data = data_long_titled[data_long_titled$GEOCODE2 %in% bgy_universe_none_maj_2016, ])

# total overlap (binary only)
mod.ipbin.titled.confirmed.total <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                         data = data_long_titled[data_long_titled$GEOCODE2 %in% bgy_universe_none_total_2016, ])

# Create table
numobs <- c(nobs(mod.ip.universe.confirmed.maj) / 2,
            nobs(mod.ipbin.universe.confirmed.maj) / 2,
            nobs(mod.ipbin.universe.confirmed.total) / 2,
            nobs(mod.ip.titled.confirmed.maj) / 2,
            nobs(mod.ipbin.titled.confirmed.maj) / 2,
            nobs(mod.ipbin.titled.confirmed.total) / 2
)
star <- stargazer(mod.ip.universe.confirmed.maj,
                  mod.ipbin.universe.confirmed.maj,
                  mod.ipbin.universe.confirmed.total,
                  mod.ip.titled.confirmed.maj,
                  mod.ipbin.titled.confirmed.maj,
                  mod.ipbin.titled.confirmed.total,
                  
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Indigenous Prop.",
                  font.size = "footnotesize",
                  column.labels = c("Eligible", "Eligible", "Eligible", "Titled", "Titled", "Titled"),
                  covariate.labels = c("Titled (Prop.)",
                                       "Titled (Binary)"),
                  add.lines = list(c("Observations", numobs),
                                   c("Excludes partially titled", rep(c("N", "N", "Y"), 2))),
                  label = "tab:idconfirmed",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling and Indigenous Identification - Source Agreement")
cat(star, sep = '\n',
    file = "id_confirmed.tex"
    )

# Table A.12: Land Titling (Continuous) and Indigenous Identification - Excluding ARMM

# Subset to exclude ARMM
data_long_rural_noarmm <- data_long[data_long$rural == 1 & data_long$reg_name != "CAR" & data_long$reg_name != "ARMM", ]
data_long_matched_noarmm <- data_long[data_long$GEOCODE2 %in% matched_samp_bgys & data_long$reg_name != "ARMM", ]
data_long_universe_noarmm <- data_long[data_long$universe_ncip_updated == 1 & data_long$reg_name != "CAR" & data_long$reg_name != "ARMM", ]
data_long_titled_noarmm <- data_long[data_long$Approved_CADT_NCIP_database == 1 & data_long$Approved_CADT_NCIP_database_munonly == 0 & data_long$reg_name != "CAR" & data_long$reg_name != "ARMM", ]

# Continuous
# All rural
mod.ip.rural.noarmm <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                     data = data_long_rural_noarmm)
# Matched
mod.ip.matched.noarmm <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                       data = data_long_matched_noarmm)
# Universe
mod.ip.universe.noarmm <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                        data = data_long_universe_noarmm)
# Titled
mod.ip.titled.noarmm <- felm(ip ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                      data = data_long_titled_noarmm)

# Binary
# All rural
mod.ipbin.rural.noarmm <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                            data = data_long_rural_noarmm)
# Matched
mod.ipbin.matched.noarmm <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                              data = data_long_matched_noarmm)
# Universe
mod.ipbin.universe.noarmm <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                               data = data_long_universe_noarmm)
# Titled
mod.ipbin.titled.noarmm <- felm(ip ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                             data = data_long_titled_noarmm)

# Create table
numobs <- c(nobs(mod.ip.rural.noarmm) / 2,
            nobs(mod.ip.universe.noarmm) / 2,
            nobs(mod.ip.matched.noarmm) / 2,
            nobs(mod.ip.titled.noarmm) / 2)
star <- stargazer(mod.ip.rural.noarmm,
                  mod.ip.universe.noarmm,
                  mod.ip.matched.noarmm,
                  mod.ip.titled.noarmm,
                  mod.ipbin.rural.noarmm,
                  mod.ipbin.universe.noarmm,
                  mod.ipbin.matched.noarmm,
                  mod.ipbin.titled.noarmm,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Indigenous Prop.",
                  font.size = "footnotesize",
                  column.labels = rep(c("All Rural", "Eligible", "Matched", "Titled"), 2),
                  covariate.labels = c("Titled (Prop.)",
                                       "Titled (Binary)"),
                  add.lines = list(c("Observations", rep(numobs, 2))),
                  label = "tab:idnoarmm",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling and Indigenous Identification - Excluding ARMM")
cat(star, sep = '\n',
    file = "id_noarmm.tex"
    )

# Table A.13: Land Titling (Binary) and Birth Registration (Two-Way FE)

# All rural
mod.bregbin.rural <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                          data = data_long4_rural)
mod.bregbin.rural2 <- felm(birth_reg ~ cadt_bin_update + have_bhs + have_street_pattern + have_highway + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                           data = data_long4_controls_rural)
# Matched
mod.bregbin.matched <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                            data = data_long4_matched)
mod.bregbin.matched2 <- felm(birth_reg ~ cadt_bin_update + have_bhs + have_street_pattern + have_highway + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                             data = data_long4_controls_matched)
# Universe
mod.bregbin.universe <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                             data = data_long4_universe)
mod.bregbin.universe2 <- felm(birth_reg ~ cadt_bin_update + have_bhs + have_street_pattern + have_highway + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                              data = data_long4_controls_universe)
# Titled
mod.bregbin.titled <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                           data = data_long4_titled)
mod.bregbin.titled2 <- felm(birth_reg ~ cadt_bin_update + have_bhs + have_street_pattern + have_highway + logpop + age + hh_size | GEOCODE2 + time | 0 | GEOCODE2,
                            data = data_long4_controls_titled)

# Create table
numobs <- c(nobs(mod.bregbin.rural) / 4,
            nobs(mod.bregbin.rural2) / 4,
            nobs(mod.bregbin.universe) / 4,
            nobs(mod.bregbin.universe2) / 4,
            nobs(mod.bregbin.matched) / 4,
            nobs(mod.bregbin.matched2) / 4,
            nobs(mod.bregbin.titled) / 4,
            nobs(mod.bregbin.titled2) / 4)
star <- stargazer(mod.bregbin.rural,
                  mod.bregbin.rural2,
                  mod.bregbin.universe,
                  mod.bregbin.universe2,
                  mod.bregbin.matched,
                  mod.bregbin.matched2,
                  mod.bregbin.titled,
                  mod.bregbin.titled2,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Birth Registration",
                  font.size = "footnotesize",
                  column.labels = c("All Rural", "All Rural", "Eligible", "Eligible", "Matched", "Matched", "Titled", "Titled"),
                  covariate.labels = c("Titled (Binary)",
                                       "Bgy. Health Center",
                                       "Street Pattern",
                                       "Highway Access",
                                       "Log Population",
                                       "Mean Age",
                                       "Mean HH Size"),
                  add.lines = list(c("Observations", numobs)),
                  label = "tab:bregpanelbin",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling (Binary) and Birth Registration (Two-Way FE)")
cat(star, sep = '\n',
    file = "breg_bin.tex"
    )

# Table A.14: Land Titling and Birth Registration (Two-Way FE) - Excluding 90% Baseline

data_long4_nohi <- reshape(bgy[-which(bgy$birth_reg_2000 > 0.90), ],
                     idvar = "GEOCODE2",
                     varying = varyingvars_panel,
                     direction = "long",
                     sep = "")
# Drop barangays that that aren't present in both census years
data_long4_nohi_nona <- data_long4_nohi[!is.na(data_long4_nohi$birth_reg), ] 
data_long4_nohi <- make.pbalanced(data_long4_nohi_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))
# Create subsets for analysis
data_long4_nohi_rural <- data_long4_nohi[which(data_long4_nohi$rural == 1 & data_long4_nohi$reg_name != "CAR"), ]
data_long4_nohi_matched <- data_long4_nohi[data_long4_nohi$GEOCODE2 %in% matched_samp_bgys, ]
data_long4_nohi_universe <- data_long4_nohi[data_long4_nohi$universe_ncip_updated == 1 & data_long4_nohi$reg_name != "CAR", ]
data_long4_nohi_titled <- data_long4_nohi[data_long4_nohi$Approved_CADT_NCIP_database == 1 & data_long4_nohi$Approved_CADT_NCIP_database_munonly == 0 & data_long4_nohi$reg_name != "CAR", ]

# Continuous
# All rural
mod.breg.nohi.rural <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                            data = data_long4_nohi_rural)
# Matched
mod.breg.nohi.matched <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                              data = data_long4_nohi_matched)
# Universe
mod.breg.nohi.universe <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                               data = data_long4_nohi_universe)
# Titled
mod.breg.nohi.titled <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                             data = data_long4_nohi_titled)

# binary
# All rural
mod.bregbin.nohi.rural <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                               data = data_long4_nohi_rural)
# Matched
mod.bregbin.nohi.matched <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                 data = data_long4_nohi_matched)
# Universe
mod.bregbin.nohi.universe <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                  data = data_long4_nohi_universe)
# Titled
mod.bregbin.nohi.titled <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                data = data_long4_nohi_titled)
# Create table
numobs <- c(nobs(mod.breg.nohi.rural) / 4,
            nobs(mod.breg.nohi.universe) / 4,
            nobs(mod.breg.nohi.matched) / 4,
            nobs(mod.breg.nohi.titled) / 4)
star <- stargazer(mod.breg.nohi.rural,
                  mod.breg.nohi.universe,
                  mod.breg.nohi.matched,
                  mod.breg.nohi.titled,
                  mod.bregbin.nohi.rural,
                  mod.bregbin.nohi.universe,
                  mod.bregbin.nohi.matched,
                  mod.bregbin.nohi.titled,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Birth Registration",
                  font.size = "footnotesize",
                  column.labels = rep(c("All Rural", "Eligible", "Matched", "Titled"), 2),
                  covariate.labels = c("Titled (Prop.)",
                                       "Titled (Binary)"
                  ),
                  add.lines = list(c("Observations", rep(numobs, 2))),
                  label = "tab:bregnohi",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling and Birth Registration (Two-Way FE) - Excluding 90\\% Baseline")
cat(star, sep = '\n', 
    file = "breg_nohi.tex"
    )
rm(data_long4_nohi, data_long4_nohi_nona, data_long4_nohi_rural,
   data_long4_nohi_matched, data_long4_nohi_universe, data_long4_nohi_titled)

# Table A.15: Land Titling (Continuous) and Lagged Birth Registration
data_long4_sort <- arrange(data_long4, GEOCODE2, time)
data_long4_sort2 <- data_long4_sort %>% group_by(GEOCODE2) %>% mutate(lag_birth_reg = lag(birth_reg))
# Drop barangays from the first year (no lag)
data_long4_lag_nona <- data_long4_sort2[!is.na(data_long4_sort2$lag_birth_reg), ] 
data_long4_lag <- make.pbalanced(data_long4_lag_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))
# Create subsets for analysis
data_long4_lag_rural <- data_long4_lag[which(data_long4_lag$rural == 1 & data_long4_lag$reg_name != "CAR"), ]
data_long4_lag_matched <- data_long4_lag[data_long4_lag$GEOCODE2 %in% matched_samp_bgys, ]
data_long4_lag_universe <- data_long4_lag[data_long4_lag$universe_ncip_updated == 1 & data_long4_lag$reg_name != "CAR", ]
data_long4_lag_titled <- data_long4_lag[data_long4_lag$Approved_CADT_NCIP_database == 1 & data_long4_lag$Approved_CADT_NCIP_database_munonly == 0 & data_long4_lag$reg_name != "CAR", ]

# All rural
mod.breg.lag.rural <- felm(lag_birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                           data = data_long4_lag_rural)
# Matched
mod.breg.lag.matched <- felm(lag_birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                             data = data_long4_lag_matched)
# Universe
mod.breg.lag.universe <- felm(lag_birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                              data = data_long4_lag_universe)
# Titled
mod.breg.lag.titled <- felm(lag_birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                            data = data_long4_lag_titled)

# Create table
numobs <- c(nobs(mod.breg.lag.rural) / 3,
            nobs(mod.breg.lag.universe) / 3,
            nobs(mod.breg.lag.matched) / 3,
            nobs(mod.breg.lag.titled) / 3)
star <- stargazer(mod.breg.lag.rural,
                  mod.breg.lag.universe,
                  mod.breg.lag.matched,
                  mod.breg.lag.titled,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Lagged Birth Registration",
                  font.size = "footnotesize",
                  column.labels = c("All Rural", "Eligible", "Matched", "Titled"),
                  covariate.labels = c("Titled Prop."),
                  add.lines = list(c("Observations", numobs)),
                  label = "tab:breglag",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling (Continuous) and Lagged Birth Registration (Two-Way FE)")
cat(star, sep = '\n',
    file = "breg_lag.tex"
    )
rm(data_long4_sort, data_long4_sort2, data_long4_lag_nona,
   data_long4_lag, data_long4_lag_rural, data_long4_lag_matched,
   data_long4_lag_universe, data_long4_lag_titled) # for memory


# Table A.16: Land Titling and Birth Registration (Quasi-RDD)

bgy$cadt_approval_year_first_update <- bgy$cadt_approval_year_first
bgy$cadt_approval_year_first_update[bgy$Approved_CADT_NCIP_database_munonly == 1] <- NA

# Change from 2000 to 2010
bgy$birth_reg_delta <- bgy$birth_reg_2010 - bgy$birth_reg_2000
# Change from 2000 to 2007
bgy$birth_reg_delta_2007 <- bgy$birth_reg_2007 - bgy$birth_reg_2000

# Distance from titling year to 2010 census
bgy$census_distance2010 <- bgy$cadt_approval_year_first_update - 2010
bgy$census_distance2007 <- bgy$cadt_approval_year_first_update - 2007

# Subset to titled only
bgy_titled <- bgy[which(bgy$Approved_CADT_NCIP_database == 1 & bgy$Approved_CADT_NCIP_database_munonly == 0 & bgy$reg_name != "CAR"), ]

# 2010 - 8-year bandwidth
mod.quasi.rdd.breg10.all8 <- lm(birth_reg_delta ~ cadt_bin_update2010*census_distance2010, data = bgy_titled[bgy_titled$census_distance2010 %in% c(-8:8), ])

# 2010 - 4-year bandwidth
mod.quasi.rdd.breg10.all4 <- lm(birth_reg_delta ~ cadt_bin_update2010*census_distance2010, data = bgy_titled[bgy_titled$census_distance2010 %in% c(-4:4), ])

# 2007 - 4-year bandwidth
mod.quasi.rdd.breg07.all4 <- lm(birth_reg_delta_2007 ~ cadt_bin_update2007*census_distance2007, data = bgy_titled[bgy_titled$census_distance2007 %in% c(-4:4), ])

# Create table
stargazer(mod.quasi.rdd.breg07.all4,
          mod.quasi.rdd.breg10.all4,
          mod.quasi.rdd.breg10.all8,
          covariate.labels = c("Titled (2007)",
                               "Census Distance (2007)",
                               "Titled x Census Dist. (2007)",
                               "Titled (2010)",
                               "Census Distance (2010)",
                               "Titled x Census Dist. (2010)",
                               "Constant"),
          omit.stat = c("ser", "f"),
          add.lines = list(c("Bandwidth", 4, 4, 8)),
          dep.var.labels = c("Prop. Registered $\\Delta$ 2000-2007", "Prop. Registered $\\Delta$ 2000-2010"),
          title = "Land Titling and Birth Registration (Quasi-RDD)",
          label = "tab:bregqrdd",
          out = "breg_rdd_bin.tex"
          )

# Table A.17: Land Titling and Birth Registration (First-difference with filing year)

# Restrict to barangays associated with a CADT not missing a filing year
bgy_filed <- bgy[!is.na(bgy$filed_year) & bgy$reg_name != "CAR", ]

# control for filed year as a continuous variable
breg.delta.filed <- lm(birth_reg_delta ~ cadt_prop2010 + filed_year, data = bgy_filed)

# filing year fixed effects
breg.delta.filed2 <- lm(birth_reg_delta ~ cadt_prop2010 + as.factor(filed_year), data = bgy_filed)

# repeat with binary titling variable
bregbin.delta.filed <- lm(birth_reg_delta ~ cadt_bin_update2010 + filed_year, data = bgy_filed)

bregbin.delta.filed2 <- lm(birth_reg_delta ~ cadt_bin_update2010 + as.factor(filed_year), data = bgy_filed)

# Create table
star <- stargazer(breg.delta.filed,
                  breg.delta.filed2,
                  bregbin.delta.filed,
                  bregbin.delta.filed2,
                  omit = c("as.factor",
                           "Constant"),
                  omit.stat = c("ser", "f"),
                  dep.var.labels = "Registered Prop. $\\Delta$ 2000-2010",
                  covariate.labels = c("Titled Prop.",
                                       "Titled (Binary)",
                                       "Filed Year"),
                  label = "tab:bregfile",
                  add.lines = list(c("Filed Year FE", "N", "Y", "N", "Y")),
                  title = "Land Titling and Birth Registration (First-difference with filing year)")
cat(star, 
    sep = "\n", 
    file = "breg_filing.tex"
    )

# Table A.18: Land Titling and Birth Registration (First-difference with Province and Region Fixed Effects)

bgy$birth_reg_delta <- bgy$birth_reg_2010 - bgy$birth_reg_2000

bgy_rural <- bgy[which(bgy$rural == 1 & bgy$reg_name != "CAR"), ]
bgy_universe <- bgy[which(bgy$universe_ncip_updated == 1 & bgy$reg_name != "CAR"), ]
bgy_matched <- bgy[which(bgy$GEOCODE2 %in% matched_samp_bgys), ]
bgy_titled <- bgy[which(bgy$Approved_CADT_NCIP_database == 1 & bgy$Approved_CADT_NCIP_database_munonly == 0 & bgy$reg_name != "CAR"), ]

# Province FE
breg.delta.prov <- lm(birth_reg_delta ~ cadt_prop2010 + I(prov_code_2000), data = bgy_rural)

breg.delta.prov.matched <- lm(birth_reg_delta ~ cadt_prop2010 + I(prov_code_2000), data = bgy_matched)

breg.delta.prov.universe <- lm(birth_reg_delta ~ cadt_prop2010 + I(prov_code_2000), data = bgy_universe)

breg.delta.prov.titled <- lm(birth_reg_delta ~ cadt_prop2010 + I(prov_code_2000), data = bgy_titled)

# Region FE
breg.delta.reg <- lm(birth_reg_delta ~ cadt_prop2010 + I(reg_name), data = bgy_rural)

breg.delta.reg.matched <- lm(birth_reg_delta ~ cadt_prop2010 + I(reg_name), data = bgy_matched)

breg.delta.reg.universe <- lm(birth_reg_delta ~ cadt_prop2010 + I(reg_name), data = bgy_universe)

breg.delta.reg.titled <- lm(birth_reg_delta ~ cadt_prop2010 + I(reg_name), data = bgy_titled)

# Create table
star <- stargazer(breg.delta.prov,
                  breg.delta.prov.universe,
                  breg.delta.prov.matched,
                  breg.delta.prov.titled,
                  breg.delta.reg,
                  breg.delta.reg.universe,
                  breg.delta.reg.matched,
                  breg.delta.reg.titled,
                  omit = c("I\\(", "Constant"),
                  omit.stat = c("rsq", "f", "ser"),
                  type = "latex",
                  font.size = "footnotesize",
                  label = "tab:bregdidprovreg",
                  column.labels = rep(c("All Rural", "Eligible", "Matched", "Titled"), 2),
                  dep.var.labels = "Birth Reg. $\\Delta$ 2000-2010",
                  covariate.labels = "Titled Prop.",
                  add.lines = list(c("Prov FE", c(rep("Y", 4), rep("N", 4))),
                                   c("Reg FE", c(rep("N", 4), rep("Y", 4)))
                  ),
                  title = "Land Titling and Birth Registration (First-difference with Province and Region Fixed Effects)")
cat(star, sep = '\n', 
    file = "breg_provreg.tex"
)

# Table A.19: Land Titling and Birth Registration (Conley SE)

mod.breg.rural.conley <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | lat + lon, data = data_long4_rural, keepCX = TRUE)
SE.breg.rural <- ConleySEs(reg = mod.breg.rural.conley,
                           unit = "GEOCODE2",
                           time = "time",
                           lat = "lat", lon = "lon",
                           dist_fn = "SH", dist_cutoff = 500,
                           lag_cutoff = 5,
                           cores = 1,
                           verbose = FALSE)
ses.breg.rural <- sapply(SE.breg.rural, function(x) diag(sqrt(x))) %>% round(3)

mod.breg.matched.conley <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | lat + lon, data = data_long4_matched, keepCX = TRUE)
SE.breg.matched <- ConleySEs(reg = mod.breg.matched.conley,
                             unit = "GEOCODE2",
                             time = "time",
                             lat = "lat", lon = "lon",
                             dist_fn = "SH", dist_cutoff = 500,
                             lag_cutoff = 5,
                             cores = 1,
                             verbose = FALSE)
ses.breg.matched <- sapply(SE.breg.matched, function(x) diag(sqrt(x))) %>% round(3)

mod.breg.universe.conley <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | lat + lon, data = data_long4_universe, keepCX = TRUE)
SE.breg.universe <- ConleySEs(reg = mod.breg.universe.conley,
                              unit = "GEOCODE2",
                              time = "time",
                              lat = "lat", lon = "lon",
                              dist_fn = "SH", dist_cutoff = 500,
                              lag_cutoff = 5,
                              cores = 1,
                              verbose = FALSE)
ses.breg.universe <- sapply(SE.breg.universe, function(x) diag(sqrt(x))) %>% round(3)

mod.breg.titled.conley <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | lat + lon, data = data_long4_titled, keepCX = TRUE)
SE.breg.titled <- ConleySEs(reg = mod.breg.titled.conley,
                            unit = "GEOCODE2",
                            time = "time",
                            lat = "lat", lon = "lon",
                            dist_fn = "SH", dist_cutoff = 500,
                            lag_cutoff = 5,
                            cores = 1,
                            verbose = FALSE)
ses.breg.titled <- sapply(SE.breg.titled, function(x) diag(sqrt(x))) %>% round(3)

ses.all.breg <- c(ses.breg.rural[3], ses.breg.universe[3], ses.breg.matched[3], ses.breg.titled[3])

# Create table
numobs <- c(nobs(mod.breg.rural.conley) / 4,
            nobs(mod.breg.universe.conley) / 4,
            nobs(mod.breg.matched.conley) / 4,
            nobs(mod.breg.titled.conley) / 4)
rsqs <- round(c(summary(mod.breg.rural.conley)$r.squared[1],
                summary(mod.breg.universe.conley)$r.squared[1],
                summary(mod.breg.matched.conley)$r.squared[1],
                summary(mod.breg.titled.conley)$r.squared[1]), 3)

star <- stargazer(mod.breg.rural.conley,
                  mod.breg.universe.conley,
                  mod.breg.matched.conley,
                  mod.breg.titled.conley,
                  se = ses.all.breg,
                  omit.stat = c("rsq", "f", "ser", "n"),
                  type = "latex",
                  font.size = "footnotesize",
                  label = "tab:bregpanelconley",
                  column.labels = c("All Rural", "Eligible", "Matched", "Titled"),
                  dep.var.labels = "Birth Registration",
                  covariate.labels = "Titled Prop.",
                  add.lines = list(c("Observations", numobs),
                                   c("$R^2$", rsqs)),
                  title = "Land Titling and Birth Registration (Conley SE)")
cat(star,
    sep = '\n',
    file = "breg_panel_conley.tex"
    )

# Table A.20: P-Values for Birth Registration Analysis (Two-Way Clustering)

# Multiway clustering wild cluster bootstrap

# function takes outcome variable, single treatment variable, dataset, number of bootstrap iterations
# returns p-value single coefficient of interest (first coefficient)
multiWCB <- function(y, tx, data, nboots){
  # estimate restricted 2FE model using feols from fixest WITH DOUBLE CLUSTERING  
  mod <- paste(y, tx, sep = "~")
  mod.feols <- feols(formula(paste(mod, " | GEOCODE2 + time")), cluster = c("GEOCODE2", "time"), data, warn = F)
  # get actual t stat
  tstar <- mod.feols$coeftable[3]
  # get actual BETA
  betastar <- mod.feols$coeftable[1]
  # get actual WALD STAT
  waldstar <- wald(mod.feols)$stat
  # fit unrestricted model
  mod.feols.ur <- feols(formula("birth_reg ~ 1 | GEOCODE2 + time"), cluster = c("GEOCODE2", "time"), data)
  # residuals from unrestricted model
  resid_ur <- mod.feols.ur$residuals
  # predicted values from unrestricted model
  pred_feols_ur <- predict(mod.feols.ur, newdata = data)
  # number of clusters on first dimension
  num_clust1 <- length(unique(data[, "GEOCODE2"]))
  # number of clusters on second dimension
  num_clust2 <- length(unique(data[, "time"]))
  # number of clusters in the intersection
  num_clust12 <- nrow(data) # only one observation in each cluster intersection
  # number of observations per cluster
  obs_per_clust1 <- nrow(data) / num_clust1
  obs_per_clust2 <- nrow(data) / num_clust2
  boot_t <- c()
  for (i in 1:nboots){
    # draw disturbance terms from a Rachemacher random variable
    # six point distribution
    sixpoint <- c(-sqrt(1.5), -sqrt(1), -sqrt(0.5), sqrt(0.5), sqrt(1), sqrt(1.5))
    draw_dist_dim1 <- sample(sixpoint, size = num_clust1, prob = rep(1/6, 6), replace = T) # first clustering dimension
    draw_dist_dim2 <- sample(sixpoint, size = num_clust2, prob = rep(1/6, 6), replace = T) # second clustering dimension
    draw_dist_dim12 <- sample(sixpoint, size = num_clust12, prob = rep(1/6, 6), replace = T) # intersection of first and second clustering dimensions     # # create vector of disturbance terms with length n
    dist_vec_dim1 <- rep(draw_dist_dim1, obs_per_clust1)
    dist_vec_dim2 <- rep(draw_dist_dim2, each = obs_per_clust2)
    dist_vec_dim12 <- draw_dist_dim12
    # multiply actual errors by disturbance terms
    new_errors_dim1 <- resid_ur * dist_vec_dim1
    new_errors_dim2 <- resid_ur * dist_vec_dim2
    new_errors_dim12 <- resid_ur * dist_vec_dim12
    # create simulated dependent variable
    sim_dv_dim1 <- pred_feols_ur + new_errors_dim1
    sim_dv_dim2 <- pred_feols_ur + new_errors_dim2
    sim_dv_dim12 <- pred_feols_ur + new_errors_dim12
    # new data frame
    new_data <- data.frame(sim_dv_dim1,
                           sim_dv_dim2,
                           sim_dv_dim12,
                           tx = data[, tx], 
                           GEOCODE2 = data$GEOCODE2, 
                           time = data$time)
    # run fixed effects regression with clustering on each dimension + the intersection
    sim.feols1 <- feols(formula("sim_dv_dim1 ~ tx | GEOCODE2 + time"), cluster = "GEOCODE2", data = new_data)
    sim.feols2 <- feols(formula("sim_dv_dim2 ~ tx | GEOCODE2 + time"), cluster = "time", data = new_data)
    sim.feols12 <- feols(formula("sim_dv_dim12 ~ tx | GEOCODE2 + time"), data = new_data)
    # get standard errors
    se.feols1 <- sim.feols1$coeftable[2]
    se.feols2 <- sim.feols2$coeftable[2]
    se.feols12 <- sim.feols12$coeftable[2]
    # calculate standard error with two way clustering
    se2w <- sqrt(se.feols1^2 + se.feols2^2 - se.feols12^3)
    # get t values
    # randomly select an estimate of beta from one of the three regressions
    betas <- c(sim.feols1$coeftable[1],
               sim.feols2$coeftable[1],
               sim.feols12$coeftable[1]
               )
    beta_select <- sample(betas, 1, prob = rep(1 / length(betas), length(betas)))
    boot_t[i] <- beta_select / se2w
  }
  # get wild bootstrap p-value
  boot_pval_t <- mean(abs(unlist(boot_t)) > abs(as.numeric(tstar)))
  boot_pval_wald <- mean(abs(unlist(boot_t)) > abs(as.numeric(waldstar)))
  return(list(boot_pval_t, boot_pval_wald))
}

set.seed(02139)
p.mwc.breg.rural <- multiWCB(y = "birth_reg", tx = "cadt_prop", data = data_long4_rural, nboots = 1000)
p.mwc.breg.matched <- multiWCB(y = "birth_reg", tx = "cadt_prop", data = data_long4_matched, nboots = 1000)
p.mwc.breg.universe <- multiWCB(y = "birth_reg", tx = "cadt_prop", data = data_long4_universe, nboots = 1000)
p.mwc.breg.titled <- multiWCB(y = "birth_reg", tx = "cadt_prop", data = data_long4_titled, nboots = 1000)

# create table comparing bootstrapped two-way clusters and one-way clusters
bootstrapped <- c(p.mwc.breg.rural[[1]],
                  p.mwc.breg.universe[[1]],
                  p.mwc.breg.matched[[1]],
                  p.mwc.breg.titled[[1]]
)
oneway <- c(summary(mod.breg.rural)$coefficients[4],
            summary(mod.breg.universe)$coefficients[4],
            summary(mod.breg.matched)$coefficients[4],
            summary(mod.breg.titled)$coefficients[4])

pcompare <- data.frame(bootstrapped, oneway)
rownames(pcompare) <- c("All Rural", "Eligible", "Matched", "Titled")
colnames(pcompare) <- c("Two-Way Clustering (Bootstrapped)", "One-Way Clustering")
pcomparex <- xtable(pcompare,
                    digits = 4,
                    caption = "P-Values for Birth Registration Analysis",
                    label = "tab:mwcbreg")
print(pcomparex,
      type = 'latex',
      file = "mwc_breg.tex"
      )

# Table A.21: Land Titling and Birth Registration - Source Agreement

# subsetting only to undisputed barangays (repeated above)
# map says over 50% overlap and list says 1
confirmed_maj_2016 <- bgy_universe$GEOCODE2[bgy_universe$cadt_prop2016 > 0.5 & bgy_universe$cadt_bin_update2016 == 1]

# map says total overlap and list says 1
confirmed_total_2016 <- bgy_universe$GEOCODE2[bgy_universe$cadt_prop2016 == 1 & bgy_universe$cadt_bin_update2016 == 1]

# map says zero overlap and list says zero overlap (in 2016)
confirmed_none_2016 <- bgy_universe$GEOCODE2[bgy_universe$cadt_prop2016 == 0 & bgy_universe$cadt_bin_update2016 == 0]

bgy_universe_none_total_2016 <- c(confirmed_total_2016, confirmed_none_2016)
bgy_universe_none_maj_2016 <- c(confirmed_maj_2016, confirmed_none_2016)

# Universe
# majority overlap - continuous
mod.breg.universe.confirmed.maj <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                                        data = data_long4_universe[data_long4_universe$GEOCODE2 %in% bgy_universe_none_maj_2016, ])
# majority overlap - binary
mod.bregbin.universe.confirmed.maj <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                           data = data_long4_universe[data_long4_universe$GEOCODE2 %in% bgy_universe_none_maj_2016, ])
# total overlap (binary only)
mod.bregbin.universe.confirmed.total <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                             data = data_long4_universe[data_long4_universe$GEOCODE2 %in% bgy_universe_none_total_2016, ])

# Titled
# majority overlap - continuous
mod.breg.titled.confirmed.maj <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                                      data = data_long4_titled[data_long4_titled$GEOCODE2 %in% bgy_universe_none_maj_2016, ])
# majority overlap - binary
mod.bregbin.titled.confirmed.maj <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                         data = data_long4_titled[data_long4_titled$GEOCODE2 %in% bgy_universe_none_maj_2016, ])
# total overlap (binary only)
mod.bregbin.titled.confirmed.total <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                           data = data_long4_titled[data_long4_titled$GEOCODE2 %in% bgy_universe_none_total_2016, ])

# Create table
numobs <- c(nobs(mod.breg.universe.confirmed.maj) / 4,
            nobs(mod.bregbin.universe.confirmed.maj) / 4,
            nobs(mod.bregbin.universe.confirmed.total) / 4,
            nobs(mod.breg.titled.confirmed.maj) / 4,
            nobs(mod.bregbin.titled.confirmed.maj) / 4,
            nobs(mod.bregbin.titled.confirmed.total) / 4
)
star <- stargazer(mod.breg.universe.confirmed.maj,
                  mod.bregbin.universe.confirmed.maj,
                  mod.bregbin.universe.confirmed.total,
                  mod.breg.titled.confirmed.maj,
                  mod.bregbin.titled.confirmed.maj,
                  mod.bregbin.titled.confirmed.total,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Birth Registration",
                  font.size = "footnotesize",
                  column.labels = c("Eligible", "Eligible", "Eligible", "Titled", "Titled", "Titled"),
                  covariate.labels = c("Titled (Prop.)",
                                       "Titled (Binary)"),
                  add.lines = list(c("Observations", numobs),
                                   c("Excludes partially titled", rep(c("N", "N", "Y"), 2))),
                  label = "tab:bregconfirmed",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling and Birth Registration - Source Agreement")
cat(star, sep = '\n', 
    file = "breg_confirmed.tex"
)

# Table A.22:  Land Titling and Birth Registration - Excluding ARMM

# Subset to exclude ARMM
data_long4_rural_noarmm <- data_long4[which(data_long4$rural == 1 & data_long4$reg_name != "CAR" & data_long4$reg_name != "ARMM"), ]
data_long4_matched_noarmm <- data_long4[data_long4$GEOCODE2 %in% matched_samp_bgys & data_long4$reg_name != "ARMM", ]
data_long4_universe_noarmm <- data_long4[data_long4$universe_ncip_updated == 1 & data_long4$reg_name != "CAR" & data_long4$reg_name != "ARMM", ]
data_long4_titled_noarmm <- data_long4[data_long4$Approved_CADT_NCIP_database == 1 & data_long4$Approved_CADT_NCIP_database_munonly == 0 & data_long4$reg_name != "CAR" & data_long4$reg_name != "ARMM", ]

# continuous
# All rural
mod.breg.rural.noarmm <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                              data = data_long4_rural_noarmm)
# Matched
mod.breg.matched.noarmm <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                                data = data_long4_matched_noarmm)
# Universe
mod.breg.universe.noarmm <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                                 data = data_long4_universe_noarmm)
# Titled
mod.breg.titled.noarmm <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                               data = data_long4_titled_noarmm)

# binary
# All rural
mod.bregbin.rural.noarmm <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                 data = data_long4_rural_noarmm)
# Matched
mod.bregbin.matched.noarmm <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                   data = data_long4_matched_noarmm)
# Universe
mod.bregbin.universe.noarmm <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                    data = data_long4_universe_noarmm)
# Titled
mod.bregbin.titled.noarmm <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                  data = data_long4_titled_noarmm)

# Create table
numobs <- c(nobs(mod.breg.rural.noarmm) / 4,
            nobs(mod.breg.universe.noarmm) / 4,
            nobs(mod.breg.matched.noarmm) / 4,
            nobs(mod.breg.titled.noarmm) / 4)
star <- stargazer(mod.breg.rural.noarmm,
                  mod.breg.universe.noarmm,
                  mod.breg.matched.noarmm,
                  mod.breg.titled.noarmm,
                  mod.bregbin.rural.noarmm,
                  mod.bregbin.universe.noarmm,
                  mod.bregbin.matched.noarmm,
                  mod.bregbin.titled.noarmm,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Birth Registration",
                  font.size = "footnotesize",
                  column.labels = rep(c("All Rural", "Eligible", "Matched", "Titled"), 2),
                  covariate.labels = c("Titled (Prop.)",
                                       "Titled (Binary)"),
                  add.lines = list(c("Observations", rep(numobs, 2))),
                  label = "tab:bregnoarmm",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling and Birth Registration - Excluding ARMM")
cat(star, sep = '\n', 
    file = "breg_noarmm.tex"
    )

# Table A.23: Pre-Treatment Birth Registration Rates by Age

# create long dataset with bgy-birth year as the unit of analysis
include_cols <- c(colnames(bgy)[grep("birth_reg\\d{4}_2000", colnames(bgy))],
                  "cadt_prop2010",
                  "cadt_bin_update2010",
                  "GEOCODE2",
                  "reg_name",
                  "rural",
                  "universe_ncip_updated",
                  "Approved_CADT_NCIP_database",
                  "Approved_CADT_NCIP_database_munonly")
# subset to 2000 data only
bgy2000 <- bgy[, which(colnames(bgy) %in% include_cols)]
colnames(bgy2000) <- gsub("_2000", "", colnames(bgy2000))

# reshape 
varyingvars_pseudo <- c(paste0("birth_reg", c(1940:2000)))

data_long_age <- reshape(bgy2000,
                     idvar = "GEOCODE2",
                     varying = paste0("birth_reg", c(1940:2000)),
                     direction = "long",
                     sep = "")

data_long_age_nona <- data_long_age[!is.na(data_long_age$birth_reg), ] 
data_long_age <- make.pbalanced(data_long_age_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))
data_long_age$age <- 2000 - data_long_age$time

# Create subsets for analysis
data_long_age_matched <- data_long_age[data_long_age$GEOCODE2 %in% matched_samp_bgys, ]
data_long_age_universe <- data_long_age[data_long_age$universe_ncip_updated == 1 & data_long_age$reg_name != "CAR", ]
data_long_age_titled <- data_long_age[data_long_age$Approved_CADT_NCIP_database == 1 & data_long_age$Approved_CADT_NCIP_database_munonly == 0 & data_long_age$reg_name != "CAR", ]

# Regress birth registration on age interacted with 2010 titled dummy
mod.pretrend.breg.matched <- lm(birth_reg ~ age*cadt_bin_update2010, data = data_long_age_matched)
mod.pretrend.breg.matched.c <- coeftest(mod.pretrend.breg.matched, vcov. = cluster.vcov(mod.pretrend.breg.matched, cluster = data_long_age_matched$GEOCODE2))

mod.pretrend.breg.universe <- lm(birth_reg ~ age*cadt_bin_update2010, data = data_long_age_universe)
mod.pretrend.breg.universe.c <- coeftest(mod.pretrend.breg.universe, vcov. = cluster.vcov(mod.pretrend.breg.universe, cluster = data_long_age_universe$GEOCODE2))

mod.pretrend.breg.titled <- lm(birth_reg ~ age*cadt_bin_update2010, data = data_long_age_titled)
mod.pretrend.breg.titled.c <- coeftest(mod.pretrend.breg.titled, vcov. = cluster.vcov(mod.pretrend.breg.titled, cluster = data_long_age_titled$GEOCODE2))

# Create table
# get number of unique barangays for each model
mod.mat.matched <- data.frame(model.matrix(mod.pretrend.breg.matched))
mod.mat.matched$GEOCODE2 <- gsub("\\.\\d{4}", "", rownames(mod.mat.matched))

mod.mat.universe <- data.frame(model.matrix(mod.pretrend.breg.universe))
mod.mat.universe$GEOCODE2 <- gsub("\\.\\d{4}", "", rownames(mod.mat.universe))

mod.mat.titled <- data.frame(model.matrix(mod.pretrend.breg.titled))
mod.mat.titled$GEOCODE2 <- gsub("\\.\\d{4}", "", rownames(mod.mat.titled))

numobs <- c(length(unique(mod.mat.universe$GEOCODE2)),
            length(unique(mod.mat.matched$GEOCODE2)),
            length(unique(mod.mat.titled$GEOCODE2)))

star <- stargazer(mod.pretrend.breg.universe.c,
                  mod.pretrend.breg.matched.c,
                  mod.pretrend.breg.titled.c,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Birth Registration in 2000",
                  column.labels = c("Eligible", "Matched", "Titled"),
                  covariate.labels = c(
                    "Age",
                    "Titled 2010 (Binary)",
                    "Age x Titled",
                    "Constant"),
                  add.lines = list(c("Observations (barangays)", numobs)),
                  label = "tab:bregpretrends",
                  notes.label = "CRSE at bgy level",
                  title = "Pre-Treatment Birth Registration Rates by Age")
cat(star, sep = '\n',
    file = "breg_pretrends.tex"
    )
rm(bgy2000, data_long_age, data_long_age_nona, data_long_age_matched,
   data_long_age_universe, data_long_age_titled) # for memory

# Table A.24: Land Titling and 1990-2000 Demographic Trends (Eligible Universe)

# create pre-trend vars
bgy$leg_pretrend <- bgy$leg_whipple2000 - bgy$leg_whipple1990
bgy$logpop_pretrend <- bgy$logpop2000 - bgy$logpop1990
bgy$own_lot_pretrend <- bgy$own_lot2000 - bgy$own_lot1990
bgy$hs_grad_pretrend <- bgy$hs_grad2000 - bgy$hs_grad1990
bgy$strong_mat_pretrend <- bgy$strong_mat_index2000 - bgy$strong_mat_index1990

# Subsets
bgy_rural <- bgy[which(bgy$rural == 1 & bgy$reg_name != "CAR"), ]
bgy_universe <- bgy[which(bgy$universe_ncip_updated == 1 & bgy$reg_name != "CAR"), ]
bgy_matched <- bgy[which(bgy$GEOCODE2 %in% matched_samp_bgys), ]
bgy_titled <- bgy[which(bgy$Approved_CADT_NCIP_database == 1 & bgy$Approved_CADT_NCIP_database_munonly == 0 & bgy$reg_name != "CAR"), ]

# Legibility pre-trend
mod.legpre.universe <- lm(leg_pretrend ~ cadt_bin_update2010, data = bgy_universe)
# Population pre-trend
mod.logpoppre.universe <- lm(logpop_pretrend ~ cadt_bin_update2010, data = bgy_universe)
# Land ownership pre-trend
mod.own_lotpre.universe <- lm(own_lot_pretrend ~ cadt_bin_update2010, data = bgy_universe)
# HS graduation rates pre-trend
mod.hs_gradpre.universe <- lm(hs_grad_pretrend ~ cadt_bin_update2010, data = bgy_universe)
# Housing quality pre-trend
mod.strong_matpre.universe <- lm(strong_mat_pretrend ~ cadt_bin_update2010, data = bgy_universe)

# Create table
star <- stargazer(mod.legpre.universe,
                  mod.logpoppre.universe,
                  mod.own_lotpre.universe,
                  mod.hs_gradpre.universe,
                  mod.strong_matpre.universe,
                  type = "latex",
                  dep.var.labels = c("Legibility",
                                     "Logged Pop.",
                                     "Land Tenure",
                                     "HS Grad",
                                     "Housing Quality"),
                  covariate.labels = c("Titled in 2010"),
                  omit.stat = c("ser", "f"),
                  label = "tab:pretrendsuniverse",
                  title = "Land Titling and 1990-2000 Demographic Trends (Eligible Universe)"
)
cat(star, 
    sep = '\n', 
    file = "pretrends_universe.tex"
    )

# Table A.25: Land Titling and 1990-2000 Demographic Trends (Matched)

# Legibility pre-trend
mod.legpre.matched <- lm(leg_pretrend ~ cadt_bin_update2010, data = bgy_matched)
# Population pre-trend
mod.logpoppre.matched <- lm(logpop_pretrend ~ cadt_bin_update2010, data = bgy_matched)
# Land ownership pre-trend
mod.own_lotpre.matched <- lm(own_lot_pretrend ~ cadt_bin_update2010, data = bgy_matched)
# HS graduation rates pre-trend
mod.hs_gradpre.matched <- lm(hs_grad_pretrend ~ cadt_bin_update2010, data = bgy_matched)
# Housing quality pre-trend
mod.strong_matpre.matched <- lm(strong_mat_pretrend ~ cadt_bin_update2010, data = bgy_matched)

# Create table
star <- stargazer(mod.legpre.matched,
                  mod.logpoppre.matched,
                  mod.own_lotpre.matched,
                  mod.hs_gradpre.matched,
                  mod.strong_matpre.matched,
                  type = "latex",
                  dep.var.labels = c("Legibility",
                                     "Logged Pop.",
                                     "Land Tenure",
                                     "HS Grad",
                                     "Housing Quality"),
                  covariate.labels = c("Titled in 2010"),
                  omit.stat = c("ser", "f"),
                  label = "tab:pretrendsmatched",
                  title = "Land Titling and 1990-2000 Demographic Trends (Matched)"
)
cat(star, 
    sep = '\n', 
    file = "pretrends_matched.tex"
)

# Table A.26: Land Titling and 1990-2000 Demographic Trends (Titled)

# Legibility pre-trend
mod.legpre.titled <- lm(leg_pretrend ~ cadt_bin_update2010, data = bgy_titled)
# Population pre-trend
mod.logpoppre.titled <- lm(logpop_pretrend ~ cadt_bin_update2010, data = bgy_titled)
# Land ownership pre-trend
mod.own_lotpre.titled <- lm(own_lot_pretrend ~ cadt_bin_update2010, data = bgy_titled)
# HS graduation rates pre-trend
mod.hs_gradpre.titled <- lm(hs_grad_pretrend ~ cadt_bin_update2010, data = bgy_titled)
# Housing quality pre-trend
mod.strong_matpre.titled <- lm(strong_mat_pretrend ~ cadt_bin_update2010, data = bgy_titled)

# Create table
star <- stargazer(mod.legpre.titled,
                  mod.logpoppre.titled,
                  mod.own_lotpre.titled,
                  mod.hs_gradpre.titled,
                  mod.strong_matpre.titled,
                  type = "latex",
                  dep.var.labels = c("Legibility",
                                     "Logged Pop.",
                                     "Land Tenure",
                                     "HS Grad",
                                     "Housing Quality"),
                  covariate.labels = c("Titled in 2010"),
                  omit.stat = c("ser", "f"),
                  label = "tab:pretrendstitled",
                  title = "Land Titling and 1990-2000 Demographic Trends (Titled)"
)
cat(star, 
    sep = '\n', 
    file = "pretrends_titled.tex"
)

# Table A.27: Land Titling (Continuous) and Indigenous Identification - Land Value Interactions

data_long_land_nona <- data_long[!is.na(data_long$soil_index) & !is.na(data_long$mindep_bin) & !is.na(data_long$ip), ] 
data_long_land <- make.pbalanced(data_long_land_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))
# Create subsets for analysis
data_long_land_rural <- data_long_land[data_long_land$rural == 1 & data_long_land$reg_name != "CAR", ]
data_long_land_matched <- data_long_land[data_long_land$GEOCODE2 %in% matched_samp_bgys, ]
data_long_land_universe <- data_long_land[data_long_land$universe_ncip_updated == 1 & data_long_land$reg_name != "CAR", ]
data_long_land_titled <- data_long_land[data_long_land$Approved_CADT_NCIP_database == 1 & data_long_land$Approved_CADT_NCIP_database_munonly == 0 & data_long_land$reg_name != "CAR", ]

# Mineral deposits
mod.ip.mindep.rural <- felm(ip ~ cadt_prop + cadt_prop:mindep_bin | GEOCODE2 + time | 0 | GEOCODE2,
                            data = data_long_land_rural)

mod.ip.mindep.universe <- felm(ip ~ cadt_prop + cadt_prop:mindep_bin | GEOCODE2 + time | 0 | GEOCODE2,
                               data = data_long_land_universe)

mod.ip.mindep.matched <- felm(ip ~ cadt_prop + cadt_prop:mindep_bin | GEOCODE2 + time | 0 | GEOCODE2,
                              data = data_long_land_matched)

mod.ip.mindep.titled <- felm(ip ~ cadt_prop + cadt_prop:mindep_bin | GEOCODE2 + time | 0 | GEOCODE2,
                             data = data_long_land_titled)

# Soil quality
mod.ip.soil.rural <- felm(ip ~ cadt_prop + cadt_prop:soil_index | GEOCODE2 + time | 0 | GEOCODE2,
                          data = data_long_land_rural)

mod.ip.soil.universe <- felm(ip ~ cadt_prop + cadt_prop:soil_index | GEOCODE2 + time | 0 | GEOCODE2,
                             data = data_long_land_universe)

mod.ip.soil.matched <- felm(ip ~ cadt_prop + cadt_prop:soil_index | GEOCODE2 + time | 0 | GEOCODE2,
                            data = data_long_land_matched)

mod.ip.soil.titled <- felm(ip ~ cadt_prop + cadt_prop:soil_index | GEOCODE2 + time | 0 | GEOCODE2,
                           data = data_long_land_titled)

# Create tables
numobs <- c(nobs(mod.ip.mindep.rural) / 2,
            nobs(mod.ip.mindep.universe) / 2,
            nobs(mod.ip.mindep.matched) / 2,
            nobs(mod.ip.mindep.titled) / 2,
            nobs(mod.ip.soil.rural) / 2,
            nobs(mod.ip.soil.universe) / 2,
            nobs(mod.ip.soil.matched) / 2,
            nobs(mod.ip.soil.titled) / 2)

star <- stargazer(mod.ip.mindep.rural,
                  mod.ip.mindep.universe,
                  mod.ip.mindep.matched,
                  mod.ip.mindep.titled,
                  mod.ip.soil.rural,
                  mod.ip.soil.universe,
                  mod.ip.soil.matched,
                  mod.ip.soil.titled,
                  type = "latex",
                  dep.var.labels = "Indigenous Prop.",
                  font.size = "footnotesize",
                  column.labels = rep(c("All Rural", "Eligible", "Matched", "Titled"), 2),
                  order = c(2, 3, 1),
                  covariate.labels = c("Titled Prop. x Mineral Dep",
                                       "Titled Prop. x Soil Quality",
                                       "Titled Prop."),
                  add.lines = list(c("Observations", numobs)),
                  omit.stat = c("n", "ser", "f"),
                  label = "tab:iddidlandval",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling (Continuous) and Indigenous Identification - Land Value Interactions")
cat(star, sep = '\n', 
   file = "id_did_landvals.tex"
    )
rm(data_long_land_nona, data_long_land, data_long_land_rural,
   data_long_land_matched, data_long_land_universe, data_long_land_titled)

# Table A.28: Land Titling and Birth Registration by Birth Year (Pseudo-panel)

# Create pseudo panel for birth registration (by age)
include_cols <- c(colnames(bgy)[grep("birth_reg\\d{4}_2015", colnames(bgy))],
                  colnames(bgy)[grep("cadt_prop", colnames(bgy))],
                  colnames(bgy)[grep("cadt_bin_update", colnames(bgy))],
                  "GEOCODE2",
                  "reg_name",
                  "rural",
                  "universe_ncip_updated",
                  "Approved_CADT_NCIP_database",
                  "Approved_CADT_NCIP_database_munonly")

# for 2015
bgy2015 <- bgy[, which(colnames(bgy) %in% include_cols)]
colnames(bgy2015) <- gsub("_2015", "", colnames(bgy2015))

varyingvars_pseudo <- c(paste0("birth_reg", c(2000:2015)),
                        paste0("cadt_prop", c(2000:2015)),
                        paste0("cadt_bin_update", c(2000:2015)))

data_long_age <- reshape(bgy2015,
                     idvar = "GEOCODE2",
                     varying = varyingvars_pseudo,
                     direction = "long",
                     sep = "")
data_long_age_nona <- data_long_age[!is.na(data_long_age$birth_reg), ] 
data_long_age <- make.pbalanced(data_long_age_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))
# Create subsets for analysis
data_long_age_rural <- data_long_age[which(data_long_age$rural == 1 & data_long_age$reg_name != "CAR"), ]
data_long_age_matched <- data_long_age[data_long_age$GEOCODE2 %in% matched_samp_bgys, ]
data_long_age_universe <- data_long_age[data_long_age$universe_ncip_updated == 1 & data_long_age$reg_name != "CAR", ]
data_long_age_titled <- data_long_age[data_long_age$Approved_CADT_NCIP_database == 1 & data_long_age$Approved_CADT_NCIP_database_munonly == 0 & data_long_age$reg_name != "CAR", ]

# All rural
mod.breg.agepanel.rural <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                                data = data_long_age_rural)
# Matched
mod.breg.agepanel.matched <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                                  data = data_long_age_matched)
# Universe
mod.breg.agepanel.universe <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                                   data = data_long_age_universe)
# Titled
mod.breg.agepanel.titled <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                                 data = data_long_age_titled)

# Create table
numobs <- c(nobs(mod.breg.agepanel.rural) / 16,
            nobs(mod.breg.agepanel.universe) / 16,
            nobs(mod.breg.agepanel.matched) / 16,
            nobs(mod.breg.agepanel.titled) / 16)
star <- stargazer(mod.breg.agepanel.rural,
                  mod.breg.agepanel.universe,
                  mod.breg.agepanel.matched,
                  mod.breg.agepanel.titled,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Birth Registration by Birth Year",
                  font.size = "footnotesize",
                  column.labels = c("All Rural", "Eligible", "Matched", "Titled"),
                  covariate.labels = c("Titled (Prop.)",
                                       "Titled (Binary)"),
                  add.lines = list(c("Observations", numobs)),
                  label = "tab:bregpseudo",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling and Birth Registration by Birth Year (Pseudo-panel)")
cat(star, sep = '\n',
    file = "breg_pseudo_panel.tex"
    )
rm(bgy2015, data_long_age, data_long_age_nona, data_long_age_rural, data_long_age_matched,
   data_long_age_universe, data_long_age_titled)

# Table A.29: Land Titling and Birth Registration Among Accessible Barangays, 2007-2015

# Recreate dataset with just 2007 to 2015
varyingvars_panel_sub <- varyingvars_panel[-grep("2000", varyingvars_panel)]

data_long3 <- reshape(bgy,
                     idvar = "GEOCODE2",
                     varying = varyingvars_panel_sub,
                     direction = "long",
                     sep = "")
# Drop barangays that that aren't present in both census years
data_long3_nona <- data_long3[!is.na(data_long3$birth_reg), ] 
data_long3 <- make.pbalanced(data_long3_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))
# Subset to universe
data_long3_universe <- data_long3[data_long3$universe_ncip_updated == 1 & data_long3$reg_name != "CAR", ]
# Subset to barangays that had BOTH highway access and a health center in both 2000 and 2007
data_long3_universe_bhs_highway11 <- data_long3_universe[data_long3_universe$have_bhs_2000 == 1 & data_long3_universe$have_bhs_2007 == 1 & data_long3_universe$have_highway_2000 == 1 & data_long3_universe$have_highway_2007 == 1, ]

# Continuous
mod.breg.universe.serviced <- felm(birth_reg ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                                   data = data_long3_universe_bhs_highway11)
# Binary
mod.bregbin.universe.serviced <- felm(birth_reg ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                                      data = data_long3_universe_bhs_highway11)

# Create table
star <- stargazer(mod.breg.universe.serviced,
                  mod.bregbin.universe.serviced,
                  dep.var.labels = "Birth Registration",
                  covariate.labels = c("Titled (Prop.)",
                                       "Titled (Binary)"),
                  title = "Land Titling and Birth Registration Among Accessible Barangays, 2007-2015",
                  omit.stat = c("ser", "f"),
                  label = "tab:bregaccess",
                  type = "latex"
)
cat(star, sep = '\n', 
    file = "breg_access.tex"
    )
rm(data_long3, data_long3_nona, data_long3_universe, 
   data_long3_universe_bhs_highway11)

# Table A.30: Land Titling and NPA Influence

# Recreate dataset for years 2011-2014 (coverage of NPA data)
data_long_npa <- reshape(bgy,
                     idvar = "GEOCODE2",
                     varying = c(
                       "cadt_prop2011", "cadt_prop2012", "cadt_prop2013", "cadt_prop2014",
                       "cadt_bin_update2011", "cadt_bin_update2012", "cadt_bin_update2013", "cadt_bin_update2014",
                       "infl2011", "infl2012", "infl2013", "infl2014"),
                     direction = "long",
                     sep = "")
data_long_npa_nona <- data_long_npa[!is.na(data_long_npa$infl), ] 
data_long_npa <- make.pbalanced(data_long_npa_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))
# Subsets
data_long_npa_rural <- data_long_npa[which(data_long_npa$rural == 1 & data_long_npa$reg_name != "CAR"), ]
data_long_npa_matched <- data_long_npa[data_long_npa$GEOCODE2 %in% matched_samp_bgys, ]
data_long_npa_universe <- data_long_npa[data_long_npa$universe_ncip_updated == 1 & data_long_npa$reg_name != "CAR", ]
data_long_npa_titled <- data_long_npa[data_long_npa$Approved_CADT_NCIP_database == 1 & data_long_npa$Approved_CADT_NCIP_database_munonly == 0 & data_long_npa$reg_name != "CAR", ]

# continuous
# All rural
mod.npa.rural <- felm(infl ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                      data = data_long_npa_rural)
# Matched
mod.npa.matched <- felm(infl ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                        data = data_long_npa_matched)
# Universe
mod.npa.universe <- felm(infl ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                         data = data_long_npa_universe)
# Titled
mod.npa.titled <- felm(infl ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                       data = data_long_npa_titled)

# binary 
# All rural
mod.npabin.rural <- felm(infl ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                         data = data_long_npa_rural)
# Matched
mod.npabin.matched <- felm(infl ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                           data = data_long_npa_matched)
# Universe
mod.npabin.universe <- felm(infl ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                            data = data_long_npa_universe)
# Titled
mod.npabin.titled <- felm(infl ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                          data = data_long_npa_titled)

# Create table
numobs <- c(nobs(mod.npa.rural) / 4,
            nobs(mod.npa.universe) / 4,
            nobs(mod.npa.matched) / 4,
            nobs(mod.npa.titled) / 4,
            nobs(mod.npabin.rural) / 4,
            nobs(mod.npabin.universe) / 4,
            nobs(mod.npabin.matched) / 4,
            nobs(mod.npabin.titled) / 4)
star <- stargazer(mod.npa.rural,
                  mod.npa.universe,
                  mod.npa.matched,
                  mod.npa.titled,
                  mod.npabin.rural,
                  mod.npabin.universe,
                  mod.npabin.matched,
                  mod.npabin.titled,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "NPA Influence",
                  font.size = "footnotesize",
                  column.labels = rep(c("All Rural", "Eligible", "Matched", "Titled"), 2),
                  covariate.labels = c("Titled (Prop.)",
                                       "Titled (Binary)"
                  ),
                  add.lines = list(c("Observations", numobs)),
                  label = "tab:npa",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling and NPA Influence")
cat(star, sep = '\n',
    file = "npa_influence.tex"
    )
rm(data_long_npa, data_long_npa_nona, data_long_npa_rural,
   data_long_npa_matched, data_long_npa_universe, data_long_npa_titled)

# Table A.31: Land Titling and Legibility (Two-Way FE)

# Create version of data with 5 periods (including 1990)
bgy$cadt_prop1990 <- 0
bgy$cadt_bin_update1990 <- 0

varyingvars_panel_indiv <- c(
  "leg_whipple1990", "leg_whipple2000", "leg_whipple2007", "leg_whipple2010", "leg_whipple2015",
  "cadt_prop1990", "cadt_prop2000", "cadt_prop2007", "cadt_prop2010", "cadt_prop2015",
  "cadt_bin_update1990", "cadt_bin_update2000", "cadt_bin_update2007", "cadt_bin_update2010", "cadt_bin_update2015"
)

data_long5 <- reshape(bgy,
                     idvar = "GEOCODE2",
                     varying = varyingvars_panel_indiv,
                     direction = "long",
                     sep = "")
# Drop barangays that that aren't present all census years
data_long5_nona <- data_long5[!is.na(data_long5$leg_whipple), ] 
data_long5 <- make.pbalanced(data_long5_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))
# Create subsets for analysis
data_long5_rural <- data_long5[which(data_long5$rural == 1 & data_long5$reg_name != "CAR"), ]
data_long5_matched <- data_long5[data_long5$GEOCODE2 %in% matched_samp_bgys, ]
data_long5_universe <- data_long5[data_long5$universe_ncip_updated == 1 & data_long5$reg_name != "CAR", ]
data_long5_titled <- data_long5[data_long5$Approved_CADT_NCIP_database == 1 & data_long5$Approved_CADT_NCIP_database_munonly == 0 & data_long5$reg_name != "CAR", ]

# continuous
# All rural
mod.leg.rural <- felm(leg_whipple ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                      data = data_long5_rural)
# Matched
mod.leg.matched <- felm(leg_whipple ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                        data = data_long5_matched)
# Universe
mod.leg.universe <- felm(leg_whipple ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                         data = data_long5_universe)
# Titled
mod.leg.titled <- felm(leg_whipple ~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2,
                       data = data_long5_titled)

# binary
# All rural
mod.legbin.rural <- felm(leg_whipple ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                         data = data_long5_rural)
# Matched
mod.legbin.matched <- felm(leg_whipple ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                           data = data_long5_matched)
# Universe
mod.legbin.universe <- felm(leg_whipple ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                            data = data_long5_universe)
# Titled
mod.legbin.titled <- felm(leg_whipple ~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2,
                          data = data_long5_titled)

# Create table
numobs <- c(nobs(mod.leg.rural) / 5,
            nobs(mod.leg.universe) / 5,
            nobs(mod.leg.matched) / 5,
            nobs(mod.leg.titled) / 5)
star <- stargazer(mod.leg.rural,
                  mod.leg.universe,
                  mod.leg.matched,
                  mod.leg.titled,
                  mod.legbin.rural,
                  mod.legbin.universe,
                  mod.legbin.matched,
                  mod.legbin.titled,
                  type = "latex",
                  omit.stat = c("n", "ser"),
                  dep.var.labels = "Legibility (Whipple Index)",
                  font.size = "footnotesize",
                  column.labels = rep(c("All Rural", "Eligible", "Matched", "Titled"), 2),
                  covariate.labels = c("Titled (Prop.)",
                                       "Titled (Binary)"),
                  add.lines = list(c("Observations", numobs, numobs)),
                  label = "tab:legpanelcontbin",
                  notes.label = "CRSE at bgy level",
                  notes.append = FALSE,
                  title = "Land Titling and Legibility (Two-Way FE)")
cat(star, sep = '\n',
    file = "leg_both.tex"
    )
rm(data_long5, data_long5_nona, data_long5_rural, data_long5_matched,
   data_long5_universe, data_long5_titled)

# Table A.32: Priming Experiment Covariate Balance

# Regress treatment on covariates
mod.balance <- lm(ipra_any ~
                    female +
                    completed_elementary +
                    completed_hs +
                    catholic +
                    evangelical +
                    born_bgy +
                    age,
                  data = all_nomat)
summary(mod.balance)
Rbeta.hat <- coef(mod.balance)[-1]
RVR <- vcovHC(mod.balance, type <- 'HC0')[-1, -1]
W_obs <- as.numeric(Rbeta.hat %*% solve(RVR, Rbeta.hat))

# compute HSK robust Wald statistic and compare to a permutation distribution
female <- all_nomat$female
completed_elementary <- all_nomat$completed_elementary
completed_hs <- all_nomat$completed_hs
catholic <- all_nomat$catholic
evangelical <- all_nomat$evangelical
born_bgy <- all_nomat$born_bgy
age <- all_nomat$age

set.seed(02139)
sims <- 1000
W_sims <- c()

for(i in 1:sims){
  Z_sim <- rbinom(nrow(all_nomat), 1, 0.5)
  fit_sim <- lm(Z_sim ~ female +
                  completed_elementary +
                  completed_hs +
                  catholic +
                  evangelical +
                  born_bgy +
                  age)
  Rbeta.hat <- coef(fit_sim)[-1]
  RVR <- vcovHC(fit_sim, type <- 'HC0')[-1, -1]
  W_sims[i] <- as.numeric(Rbeta.hat %*% solve(RVR, Rbeta.hat))
  
}

# P-value
p <- mean(W_sims >= W_obs)
p

# Create table
stargazer(mod.balance,
          covariate.labels = c("Female",
                               "Completed Elementary",
                               "Completed HS",
                               "Catholic",
                               "Evangelical",
                               "Born in Bgy",
                               "Age"
          ),
          dep.var.labels = c("IPRA Prime"),
          omit.stat = c("ser", "f"),
          font.size = "footnotesize",
          title = "Priming Experiment Covariate Balance",
          type = "latex",
          label = "tab:expbalance",
          out = "exp_balance.tex"
          
)

# Table A.33: Recognition Prime and Tribal and National Identity, by Titling Status 

# Tribe top x titling status
mod.titledx.tribe.any.cov <- lm(tribe_top ~ ipra_any*(edu_hs_dm + rel_evangelical_dm + tribe_placement_dm + cadt_bin_update2018), data = all_nomat)

# Nationality top by titling status
mod.titledx.nationality.any.cov <- lm(nationality_top ~ ipra_any*(edu_hs_dm + rel_evangelical_dm + nationality_placement_dm + cadt_bin_update2018), data = all_nomat)

# Create table
stargazer(mod.titledx.tribe.any.cov,
          mod.titledx.nationality.any.cov,
          dep.var.labels = c("Tribe Top", "Nationality Top"),
          covariate.labels = c("IPRA Prime",
                               "Titled (Binary)",
                               "IPRA Prime x Titled"),
          title = "Recognition Prime and Tribal and National Identity, by Titling Status",
          omit = c("edu", "rel", "nationality", "placement", "Constant"),
          add.lines = list(c("Covariate Adjustment", rep("Y", 2))),
          type = "latex",
          omit.stat = c("f", "ser"),
          font.size = "footnotesize",
          label = "tab:nattribetitlingcov"
          ,
          out = "nat_tribe_titling_cov.tex"
          
)

# Table A.34: Covariate Balance (Targeted vs. Actual Survey Sample)

# subset to targeted and included barangays
bgy_survey_sub <- bgy[which(bgy$included_survey == 1 | bgy$included_survey_actual == 1), ]

# Create table
covars <- c("cadt_bin_update2018",
            "ip_2000",
            "birth_reg_2000",
            "logpop_2000",
            "ethfrac_2000", 
            "hs_grad_2000",
            "area_sqkm", "elev_mean", "elev_sd",
            "slope_mean", "soil_index", "mindep_bin", "log_ncip_dist", "log_coast_dist", "log_road_dist", "religion_catholic_2000",
            "strong_mat_roof_2000", "strong_mat_outerwalls_2000", "own_lot_2000", "have_street_pattern_2000", "have_highway_2000", "have_church_2000",
            "have_market_2000", "have_elem_2000", "have_bhs_2000", "have_water_sys_2000", "leg_whipple_2000")
covars_names <- c("Titled 2018",
                  "Indigenous Prop. 2000",
                  "Birth Registration 2000",
                  "Log. Population 2000",
                  "Ethnic Frac. 2000",
                  "Graduated HS Prop. 2000",
                  "Area (sq. km)",
                  "Elevation Mean",
                  "Elevation Std. Dev.",
                  "Slope Mean",
                  "Soil Quality Index",
                  "Mineral Deposits",
                  "Log. NCIP Dist",
                  "Log. Coast Dist",
                  "Log. Road Dist",
                  "Catholic Prop. 2000",
                  "Roof Strong Materials Prop. 2000",
                  "Wall Strong Materials Prop. 2000",
                  "Own Lot Prop. 2000",
                  "Street Pattern 2000",
                  "Highway Access 2000",
                  "Church 2000",
                  "Market 2000",
                  "Elementary 2000",
                  "Bgy. Health Center 2000",
                  "Water System 2000",
                  "Legibility (Whipple) 2000")
baltab_survey <- baltest.collect(MatchBalance(formula(paste("included_survey_actual ~", paste0(covars, collapse = "+"), sep = " " )), data = bgy_survey_sub), var.names = covars, after = FALSE)
baltab_survey_sub <- baltab_survey[, c(1, 2, 6)]
colnames(baltab_survey_sub) <- c("Included Mean", "Excluded Mean", "T Pval")
rownames(baltab_survey_sub) <- covars_names

baltab_surveyx <- xtable(baltab_survey_sub, 
                         caption = "Covariate Balance  (Targeted vs. Actual Survey Sample)",
                         digits = 3,
                         fontsize = "tiny",
                         label = "tab:samplebalance")
print(baltab_surveyx,
      type = "latex",
      caption.placement = "top",
      file = "sample_balance.tex"
)

# Table A.35: Covariate Balance (Survey Barangays vs. National Eligible Universe)

bgy_universe <- bgy[bgy$universe_ncip_updated == 1 & bgy$reg_name != "CAR", ]
# repeat rows for included barangays (so the tables compare included barangays to the whole sample, including them)
bgy_survey <- bgy[bgy$included_survey_actual == 1, ]
bgy_survey$included_survey_actual <- 0
n_universe <- nrow(bgy_universe)
bgy_universe <- rbind(bgy_universe, bgy_survey)
n_survey_universe <- nrow(bgy_universe[bgy_universe$included_survey_actual == 1, ])

# survey sample vs eligible universe
baltab_survey <- baltest.collect(MatchBalance(formula(paste("included_survey_actual ~", paste0(covars, collapse = "+"), sep = " " )), data = bgy_universe), var.names = covars_names, after = FALSE)
baltab_survey_sub <- baltab_survey[, c(1, 2, 6)]
colnames(baltab_survey_sub) <- c(paste0("Survey Mean\n (N = ", n_survey_universe, ")"), paste0("National Mean\n (N = ", n_universe, ")"), "T Pval")
rownames(baltab_survey_sub) <- covars_names

baltab_surveyx <- xtable(baltab_survey_sub, 
                         caption = "Covariate Balance (Survey Barangays vs. National Eligible Universe)",
                         digits = 3,
                         fontsize = "tiny",
                         label = "tab:surveyvsuniverse")
print(baltab_surveyx,
      type = "latex",
      caption.placement = "top",
      file = "survey_vs_universe.tex"
)

# Table A.36: Covariate Balance (Titled Survey Barangays vs. National Titled Barangays)

bgy_titled <- bgy[bgy$Approved_CADT_NCIP_database == 1 & bgy$Approved_CADT_NCIP_database_munonly == 0 & bgy$reg_name != "CAR", ]
n_titled <- nrow(bgy_titled)
bgy_titled_survey <- bgy_titled[bgy_titled$included_survey_actual == 1, ] # titled barangays in the survey sample
n_survey_titled <- nrow(bgy_titled_survey)
bgy_titled_survey$included_survey_actual <- 0
bgy_titled <- rbind(bgy_titled, bgy_titled_survey)

# titled barangays in survey sample compared to titled barangays overall
baltab_survey <- baltest.collect(MatchBalance(formula(paste("included_survey_actual ~", paste0(covars[-1], collapse = "+"), sep = " " )), data = bgy_titled), var.names = covars_names[-1], after = FALSE)
baltab_survey_sub <- baltab_survey[, c(1, 2, 6)]
colnames(baltab_survey_sub) <- c(paste0("Survey Mean\n (N = ", n_survey_titled, ")"), paste0("National Mean\n (N = ", n_titled, ")"), "T Pval")
rownames(baltab_survey_sub) <- covars_names[-1]

baltab_surveyx <- xtable(baltab_survey_sub, 
                         caption = "Covariate Balance (Titled Survey Barangays vs. National Titled Barangays)",
                         digits = 3,
                         fontsize = "tiny",
                         label = "tab:surveyvstitled")
print(baltab_surveyx,
      caption.placement = "top",
      type = "latex",
      file = "survey_vs_titled.tex"
)

# Table A.37: Recognition Primes and Tribal Identity

# demeaned covariates
all$edu_hs_dm <- all$edu_hs - mean(all$edu_hs, na.rm = T)
all$rel_evangelical_dm <- all$rel_evangelical - mean(all$rel_evangelical, na.rm = T)

all$tribe_top_placement_dm <- all$tribe_top_placement - mean(all$tribe_top_placement, na.rm = T)
all$nationality_top_placement_dm <- all$nationality_top_placement - mean(all$nationality_top_placement, na.rm = T)
all$religion_top_placement_dm <- all$religion_top_placement - mean(all$religion_top_placement, na.rm = T)
all$gender_top_placement_dm <- all$gender_top_placement - mean(all$gender_top_placement, na.rm = T)

all$tribe_placement_dm <- all$tribe_placement - mean(all$tribe_placement, na.rm = T)
all$nationality_placement_dm <- all$nationality_placement - mean(all$nationality_placement, na.rm = T)
all$religion_placement_dm <- all$religion_placement - mean(all$religion_placement, na.rm = T)
all$gender_placement_dm <- all$gender_placement - mean(all$gender_placement, na.rm = T)

# Pooled treatment
mod.tribe.any <- lm(tribe_top ~ ipra_any, data = all)
mod.tribe.any.cov <- lm(tribe_top ~ ipra_any*(edu_hs_dm + rel_evangelical_dm + tribe_top_placement_dm), data = all)

# Separate treatments
mod.tribe.type <- lm(tribe_top ~ ipra_factor, data = all)
mod.tribe.type.cov <- lm(tribe_top ~ ipra_factor*(edu_hs_dm + rel_evangelical_dm + tribe_top_placement_dm), data = all)

# Create table
stargazer(mod.tribe.any,
          mod.tribe.any.cov,
          mod.tribe.type,
          mod.tribe.type.cov,
          title = "Recognition Primes and Tribal Identity",
          dep.var.labels = "Tribe Top",
          covariate.labels = c("IPRA Prime (Any)", "Identity Prime", "Material Prime"),
          omit = c(":", "edu", "rel", "tribe", "Constant"),
          omit.stat = c("f", "ser"),
          add.lines = list(c("Covariate Adjustment", "N", "Y", "N", "Y")),
          type = "latex",
          label = "tab:tribemain",
          font.size = "footnotesize",
          out = "id_pap.tex"

)

# Wald test for difference in effects of the two priming treatments
V <- vcov(mod.tribe.type.cov)
se <- sqrt(V[2, 2] + V[3, 3] - 2*V[2, 3])
tstat <- (mod.tribe.type.cov$coefficients[2] - mod.tribe.type.cov$coefficients[3]) / se
pval <- 2*(1- pt(abs(tstat), mod.tribe.type.cov$df.residual))
pval

# Table A.38: Candidate Conjoint Results by Experimental Prime

# Pooled treatment
mod.choice.any <- lm(choice ~ (gift + promises + mayor + leader + commwork)*prime_any, data = cjoint)
mod.choice.any.c <- coeftest(mod.choice.any, vcov = cluster.vcov(mod.choice.any, cjoint$respondent_id))

# Separate treatments
mod.choice.both <- lm(choice ~ (gift + promises + mayor + leader + commwork)*ipra_factor, data = cjoint)
mod.choice.both.c <- coeftest(mod.choice.both, vcov = cluster.vcov(mod.choice.both, cjoint$respondent_id))

# Create table
covlabs <- c("Gift",
             "PG Promises",
             "Mayor Alignment",
             "Leader Endorsement",
             "Community Work",
             "IPRA Prime (Any)",
             "Gift x IPRA",
             "Promises x IPRA",
             "Mayor x IPRA",
             "Leader x IPRA",
             "Comm x IPRA",
             "Identity Prime",
             "Material Prime",
             "Gift x Identity",
             "Gift x Naterial",
             "Promises x Identity",
             "Promises x Material",
             "Mayor x Identity",
             "Mayor x Material",
             "Leader x Identity",
             "Leader x Material",
             "Comm x Identity",
             "Comm x Material",
             "Constant"
)

nobs <- c(nobs(mod.choice.any),
          nobs(mod.choice.both))
nresp <- rep(length(unique(cjoint$respondent_id)), 2)

stargazer(mod.choice.any.c,
          mod.choice.both.c,
          dep.var.labels = "Choice",
          covariate.labels = covlabs,
          no.space = T,
          title = "Candidate Conjoint Results by Experimental Prime",
          font.size = "footnotesize",
          omit = c("as.factor"),
          label = "tab:cjointpap",
          type = "latex",
          add.lines = list(c("N Obs", nobs),
                             c("N Respondents", nresp)),
          out = "cjoint_pap.tex"

)

# Wald test for difference between the effects of the two primes
V <- cluster.vcov(mod.choice.both, cjoint$respondent_id)
se <- sqrt(V[15, 15] + V[16, 16] - 2*V[15, 16]) # should be no covariance because of random assignment
tstat <- (mod.choice.both$coefficients[15] - mod.choice.both$coefficients[16]) / se
pval <- (1 - pt(abs(tstat), mod.choice.both$df.residual))*2
pval

# Table A.39: Priming Treatment and Top Ranking of Additional Identity Attributes

# Nationality
mod.nationality.any <- lm(nationality_top ~ ipra_any, data = all)
mod.nationality.any.cov <- lm(nationality_top ~ ipra_any*(edu_hs_dm + rel_evangelical_dm + nationality_top_placement_dm), data = all)

mod.nationality.type <- lm(nationality_top ~ ipra_factor, data = all)
mod.nationality.type.cov <- lm(nationality_top ~ ipra_factor*(edu_hs_dm + rel_evangelical_dm + nationality_top_placement_dm), data = all)

# Gender
mod.gender.any <- lm(gender_top ~ ipra_any, data = all)
mod.gender.any.cov <- lm(gender_top ~ ipra_any*(edu_hs_dm + rel_evangelical_dm + gender_top_placement_dm), data = all)

mod.gender.type <- lm(gender_top ~ ipra_factor, data = all)
mod.gender.type.cov <- lm(gender_top ~ ipra_factor*(edu_hs_dm + rel_evangelical_dm + gender_top_placement_dm), data = all)

# Religion
mod.religion.any <- lm(religion_top ~ ipra_any, data = all)
mod.religion.any.cov <- lm(religion_top ~ ipra_any*(edu_hs_dm + rel_evangelical_dm + religion_top_placement_dm), data = all)

mod.religion.type <- lm(religion_top ~ ipra_factor, data = all)
mod.religion.type.cov <- lm(religion_top ~ ipra_factor*(edu_hs_dm + rel_evangelical_dm + religion_top_placement_dm), data = all)

# Create table
stargazer(
  mod.nationality.any,
  mod.nationality.any.cov,
  mod.nationality.type,
  mod.nationality.type.cov,
  mod.gender.any,
  mod.gender.any.cov,
  mod.gender.type,
  mod.gender.type.cov,
  mod.religion.any,
  mod.religion.any.cov,
  mod.religion.type,
  mod.religion.type.cov,
  font.size = "footnotesize",
  title = "Priming Treatment and Top Ranking of Additional Identity Attributes",
  omit = c(":", "edu", "rel", "nationality", "gender", "Constant"),
  omit.stat = c("f", "ser"),
  covariate.labels = c("IPRA Prime (Any)", "Identity Prime", "Material Prime"),
  dep.var.labels = c("Nationality Top", "Gender Top", "Religion Top"),
  add.lines = list(c("Covariate Adjustment", rep(c("N", "Y"), 6))),
  type = "latex",
  label = "tab:allnotribe",
  column.sep.width = "1pt",
   out = "allnotribe.tex"
)

# Table A.42: Priming Effects on State Attitudes

# Create index of attitudes toward the state
state_vars <- c(
  "intentions_bgy",
  "intentions_mun",
  "capabilities_bgy",
  "capabilities_mun",
  "trust_bgy",
  "trust_mun")
# Function to impute overall sample means
imputeFun <- function(data, varnames){
  newdata <- data.frame(matrix(NA, ncol = length(varnames), nrow = nrow(data)))
  colnames(newdata) <- paste(varnames, "imp", sep = "_")
  for(i in 1:length(varnames)){
    mycol <- data[, varnames[i]]
    mycol_mean <- mean(mycol, na.rm = T)
    newdata[, i] <- ifelse(is.na(mycol) == TRUE, mycol_mean, mycol)
  }
  return(newdata)
}

state_imp <- imputeFun(all, state_vars)
prcomp_state <- prcomp(state_imp, scale = TRUE)

all$state_attitudes_index <- as.matrix(state_imp) %*% as.matrix(prcomp_state$rotation[, 1])*-1

# Experimental effects on attitudes toward the state
mod.prime.state_attitudes_index_any <- lm(state_attitudes_index ~ ipra_any, data = all)
mod.prime.state_attitudes_index <- lm(state_attitudes_index ~ ipra_factor, data = all)

# Create table
stargazer(mod.prime.state_attitudes_index_any,
          mod.prime.state_attitudes_index,
          title = "Priming Effects on State Attitudes",
          dep.var.labels = "State Attitudes Index",
          omit = c("Constant"),
          covariate.labels = c("IPRA Prime (Any)", "Identity Prime", "Material Prime"),
          type = "latex",
          omit.stat = c("ser", "f"),
          label = "tab:stateattitudesexp",
          out = "state_attitudes_exp_all.tex"
)


#### APPENDIX FIGURES ####

# Figure A.1: Land titling and indigenous identification, bandwidth analysis

varyingvars_did <- c("ip2000", "ip2010",
                     "birth_reg2000", "birth_reg2010",
                     "cadt_prop2000", "cadt_prop2010",
                     "cadt_bin_update2000", "cadt_bin_update2010"
)
data_long <- reshape(bgy,
                     idvar = "GEOCODE2",
                     varying = varyingvars_did,
                     direction = "long",
                     sep = ""
)
# Drop barangays that are not present in both census years (changed borders or code)
data_long_nona <- data_long[!is.na(data_long$ip), ] 
data_long <- make.pbalanced(data_long_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))

# CADT approval year
bgy$cadt_approval_year_first_update <- bgy$cadt_approval_year_first
bgy$cadt_approval_year_first_update[bgy$Approved_CADT_NCIP_database_munonly == 1] <- NA

# Function to subset to barangays titled within a given window (in years) around the 2010 census
windowBgys <- function(x){
  lower_bound <- 2010 - x
  upper_bound <- 2010 + x
  before_bgys <- bgy$GEOCODE2[which(bgy$cadt_approval_year_first_update >= lower_bound & bgy$cadt_approval_year_first_update <= 2010 & bgy$reg_name != "CAR")]
  after_bgys <- bgy$GEOCODE2[which(bgy$cadt_approval_year_first_update > 2010 & bgy$cadt_approval_year_first_update <= upper_bound & bgy$reg_name != "CAR")]
  print(length(before_bgys) + length(after_bgys))
  return(c(before_bgys, after_bgys))
}

windowAnalysis <- function(x, mod){
  data_window <- data_long[data_long$GEOCODE2 %in% windowBgys(x), ]
  mod.window <- felm(formula(paste0(mod, "| GEOCODE2 + time | 0 | GEOCODE2")),
                                 data = data_window)
  return(summary(mod.window))
}

# Create figure
pdf("id_did_bandwidth_prop_paper.pdf")
modlist_ipprop <- lapply(2:8, function(x) windowAnalysis(x, mod = "ip ~ cadt_prop"))
res_prop <- data.frame(matrix(NA, nrow = length(modlist_ipprop), ncol = 4))
colnames(res_prop) <- c("coef", "ci_lower", "ci_upper", "nobs")
for(i in 1:length(modlist_ipprop)){
  res_prop$coef[i] <- modlist_ipprop[[i]]$coefficients[1, 1]
  se <- modlist_ipprop[[i]]$coefficients[1, 2]
  res_prop$ci_lower[i] <- res_prop$coef[i] - 1.96*se
  res_prop$ci_upper[i] <- res_prop$coef[i] + 1.96*se
  res_prop$nobs[i] <- modlist_ipprop[[i]]$N / 2
}
plot(x = c(1:length(modlist_ipprop)),
     y = res_prop$coef,
     xlim = c(0.5, length(modlist_ipprop) + 0.5),
     ylim = c(-0.15, 0.2),
     pch = 20,
     ylab = "",
     xlab = "",
     xaxt = 'n',
     cex = 1.5,
     cex.main = 0.8)
title(ylab = "Effects of titling (Prop. IP)", line = 2.5)
title(xlab = "Bandwidth around 2010 census (years)", line = 2.1)
axis(side = 1, labels = c(2:8),
     at = c(1:length(modlist_ipprop)),
     cex.axis = 1.2)
segments(x0 = c(1:length(modlist_ipprop)),
         x1 = c(1:length(modlist_ipprop)),
         y0 = res_prop$ci_lower,
         y1 = res_prop$ci_upper,
         lwd = 2)
abline(h = 0, lty = 3)
mtext(at = c(0.1, 1:length(modlist_ipprop)), side = 1, text = c("Obs.", res_prop$nobs), cex = 0.7, padj = 8)
dev.off()

# Figure A.2: Land titling and birth registration, bandwidth analysis

# Right panel: Using 2010 census as the post period

pdf("breg_did_bandwidth_prop_paper.pdf")
modlist_bregprop <- lapply(2:8, function(x) windowAnalysis(x, mod = "birth_reg ~ cadt_prop"))
res_prop <- data.frame(matrix(NA, nrow = length(modlist_bregprop), ncol = 4))
colnames(res_prop) <- c("coef", "ci_lower", "ci_upper", "nobs")
for(i in 1:length(modlist_bregprop)){
  res_prop$coef[i] <- modlist_bregprop[[i]]$coefficients[1, 1]
  se <- modlist_bregprop[[i]]$coefficients[1, 2]
  res_prop$ci_lower[i] <- res_prop$coef[i] - 1.96*se
  res_prop$ci_upper[i] <- res_prop$coef[i] + 1.96*se
  res_prop$nobs[i] <- modlist_bregprop[[i]]$N / 2
}
plot(x = c(1:length(modlist_bregprop)),
     y = res_prop$coef,
     xlim = c(0.5, length(modlist_bregprop) + 0.5),
     ylim = c(-0.15, 0.2),
     pch = 20,
     ylab = "",
     xlab = "",
     xaxt = 'n',
     cex = 1.5,
     cex.main = 0.8)
title(ylab = "Effects of titling (Prop. reg)", line = 2.5)
title(xlab = "Bandwidth around 2010 census (years)", line = 2.1)
axis(side = 1, labels = c(2:8),
     at = c(1:length(modlist_bregprop)),
     cex.axis = 1.2)
segments(x0 = c(1:length(modlist_bregprop)),
         x1 = c(1:length(modlist_bregprop)),
         y0 = res_prop$ci_lower,
         y1 = res_prop$ci_upper,
         lwd = 2)
abline(h = 0, lty = 3)
mtext(at = c(0.1, 1:length(modlist_bregprop)), side = 1, text = c("Obs.", res_prop$nobs), cex = 0.7, padj = 8)
dev.off()

# Left panel: using 2007 as the post-period

# recreate long dataset using 2007 instead of 2010
varyingvars_did_2007 <- c(
                     "birth_reg2000", "birth_reg2007",
                     "cadt_prop2000", "cadt_prop2007")
data_long_2007 <- reshape(bgy,
                     idvar = "GEOCODE2",
                     varying = varyingvars_did_2007,
                     direction = "long"
                     ,
                     sep = ""
)
# Drop barangays that that aren't present in both census years
data_long_nona <- data_long_2007[!is.na(data_long_2007$birth_reg), ] 
data_long_2007 <- make.pbalanced(data_long_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))

# revise functions
windowBgys2007 <- function(x){
  lower_bound <- 2007 - x
  upper_bound <- 2007 + x
  before_bgys <- bgy$GEOCODE2[which(bgy$cadt_approval_year_first_update >= lower_bound & bgy$cadt_approval_year_first_update <= 2007 & bgy$reg_name != "CAR")]
  after_bgys <- bgy$GEOCODE2[which(bgy$cadt_approval_year_first_update > 2007 & bgy$cadt_approval_year_first_update <= upper_bound & bgy$reg_name != "CAR")]
  print(length(before_bgys) + length(after_bgys))
  return(c(before_bgys, after_bgys))
}

windowAnalysis2007 <- function(x, mod){
  data_window <- data_long_2007[data_long_2007$GEOCODE2 %in% windowBgys2007(x), ]
  mod.window <- felm(formula(paste0(mod, "| GEOCODE2 + time | 0 | GEOCODE2")),
                     data = data_window)
  return(summary(mod.window))
}

pdf("breg_did_bandwidth_prop2007_paper.pdf")
modlist_bregprop2007 <- lapply(2:5, function(x) windowAnalysis2007(x, mod = "birth_reg ~ cadt_prop"))
res_prop <- data.frame(matrix(NA, nrow = length(modlist_bregprop2007), ncol = 4))
colnames(res_prop) <- c("coef", "ci_lower", "ci_upper", "nobs")
for(i in 1:length(modlist_bregprop2007)){
  res_prop$coef[i] <- modlist_bregprop2007[[i]]$coefficients[1, 1]
  se <- modlist_bregprop2007[[i]]$coefficients[1, 2]
  res_prop$ci_lower[i] <- res_prop$coef[i] - 1.96*se
  res_prop$ci_upper[i] <- res_prop$coef[i] + 1.96*se
  res_prop$nobs[i] <- modlist_bregprop2007[[i]]$N / 2
}

plot(x = c(1:length(modlist_bregprop2007)),
     y = res_prop$coef,
     xlim = c(0.5, length(modlist_bregprop2007) + 0.5),
     ylim = c(-0.15, 0.2),
     pch = 20,
     ylab = "",
     xlab = "",
     xaxt = 'n',
     cex = 1.5,
     cex.main = 0.8)
title(ylab = "Effects of titling (Prop. reg)", line = 2.5)
title(xlab = "Bandwidth around 2007 census (years)", line = 2.1)
axis(side = 1, labels = c(2:5),
     at = c(1:length(modlist_bregprop2007)),
     cex.axis = 1.2)
segments(x0 = c(1:length(modlist_bregprop2007)),
         x1 = c(1:length(modlist_bregprop2007)),
         y0 = res_prop$ci_lower,
         y1 = res_prop$ci_upper,
         lwd = 2)
abline(h = 0, lty = 3)
mtext(at = c(0.1, 1:length(modlist_bregprop2007)), side = 1, text = c("Obs.", res_prop$nobs), cex = 0.7, padj = 8)
dev.off()

# Figure A.3: Estimated effects of titling on barangay-level public facilities

# Create public goods index
data_long4$pg_index <- rowSums(data_long4[, c("have_bhs",
                                            "have_street_pattern",
                                            "have_highway",
                                            "have_elem",
                                            "have_water_sys"
)]) / 5
# Drop barangays that that aren't present in both census years
data_long4_pg_nona <- data_long4[!is.na(data_long4$pg_index), ] 
data_long4_pg <- make.pbalanced(data_long4_pg_nona, balance.type = "shared.individuals", index = c('GEOCODE2', 'time'))
# Create subsets for analysis
data_long4_pg_rural <- data_long4_pg[which(data_long4_pg$rural == 1 & data_long4_pg$reg_name != "CAR"), ]
data_long4_pg_matched <- data_long4_pg[data_long4_pg$GEOCODE2 %in% matched_samp_bgys, ]
data_long4_pg_universe <- data_long4_pg[data_long4_pg$universe_ncip_updated == 1 & data_long4_pg$reg_name != "CAR", ]
data_long4_pg_titled <- data_long4_pg[data_long4_pg$Approved_CADT_NCIP_database == 1 & data_long4_pg$Approved_CADT_NCIP_database_munonly == 0 & data_long4_pg$reg_name != "CAR", ]

# Function to generate plots
allsubsFun_both <- function(dv, plot_title){
  mymod <- formula(paste(dv, "~ cadt_prop | GEOCODE2 + time | 0 | GEOCODE2"))
  # mod.all <- felm(mymod, data = data_long4_pg_rural)
  mod.matched <- felm(mymod, data = data_long4_pg_matched)
  mod.universe <- felm(mymod, data = data_long4_pg_universe)
  mod.titled <- felm(mymod, data = data_long4_pg_titled)
  # binary
  mymodbin <- formula(paste(dv, "~ cadt_bin_update | GEOCODE2 + time | 0 | GEOCODE2"))
  modbin.matched <- felm(mymodbin, data = data_long4_pg_matched)
  modbin.universe <- felm(mymodbin, data = data_long4_pg_universe)
  modbin.titled <- felm(mymodbin, data = data_long4_pg_titled)
  
  # plot
  modlist <- list(mod.matched, mod.universe, mod.titled)
  modlist_bin <- list(modbin.matched, modbin.universe, modbin.titled)
  
  pdf(paste0(dv, "_plot.pdf"))
  # prop
  res_prop <- data.frame(matrix(NA, nrow = length(modlist), ncol = 3))
  colnames(res_prop) <- c("coef", "ci_lower", "ci_upper")
  for(i in 1:length(modlist)){
    res_prop$coef[i] <- summary(modlist[[i]])$coefficients[1]
    se <- summary(modlist[[i]])$coefficients[2]
    res_prop$ci_lower[i] <- res_prop$coef[i] - 1.96*se
    res_prop$ci_upper[i] <- res_prop$coef[i] + 1.96*se
  }
  # bin
  res_bin <- data.frame(matrix(NA, nrow = length(modlist_bin), ncol = 3))
  colnames(res_bin) <- c("coef", "ci_lower", "ci_upper")
  for(i in 1:length(modlist_bin)){
    res_bin$coef[i] <- summary(modlist_bin[[i]])$coefficients[1]
    se <- summary(modlist_bin[[i]])$coefficients[2]
    res_bin$ci_lower[i] <- res_bin$coef[i] - 1.96*se
    res_bin$ci_upper[i] <- res_bin$coef[i] + 1.96*se
  }
  
  # prop
  plot(x = c(1:length(modlist)) - 0.1,
       y = res_prop$coef,
       xlim = c(0.5, length(modlist) + 0.5),
       ylim = c(-0.15, 0.2),
       pch = 20,
       ylab = "Effect of titling",
       xlab = "",
       xaxt = 'n',
       main = plot_title,
       cex = 1.5,
       cex.main = 1)
  axis(side = 1, labels = c("Matched", "Eligible", "Titled"),
       at = c(1:length(modlist)),
       cex.axis = 1.2)
  segments(x0 = c(1:length(modlist)) - 0.1,
           x1 = c(1:length(modlist)) - 0.1,
           y0 = res_prop$ci_lower,
           y1 = res_prop$ci_upper,
           lwd = 2)
  abline(h = 0, lty = 3)
  # bin
  points(x = c(1:length(modlist_bin)) + 0.1,
         y = res_bin$coef,
         pch = 20,
         col = "gray",
         cex = 1.5)
  segments(x0 = c(1:length(modlist_bin)) + 0.1,
           x1 = c(1:length(modlist_bin)) + 0.1,
           y0 = res_bin$ci_lower,
           y1 = res_bin$ci_upper,
           lwd = 2,
           col = "gray")
  
  legend("topright",
         c("Continuous", "Binary"),
         lty = c(1, 1),
         lwd = c(2, 2),
         col = c("black", "gray"),
         cex = 0.8)
  dev.off()
  
  return(list(modlist,
              modlist_bin)
  )
}

allsubsFun_both("pg_index", "Public Goods Index")
allsubsFun_both("have_bhs", "Barangay Health Center")
allsubsFun_both("have_street_pattern", "Street Pattern")
allsubsFun_both("have_highway", "Highway Access")
allsubsFun_both("have_elem", "Elementary School")
allsubsFun_both("have_water_sys", "Water System")

rm(data_long4_pg, data_long4_pg_nona, data_long4_pg_rural,
   data_long4_pg_universe, data_long4_pg_matched, data_long4_pg_titled)

# Figure A.6: Effects of priming on individual components of state attitudes index 

# Create index of state attitude variables using PCA
state_vars <- c(
  "intentions_bgy",
  "intentions_mun",
  "capabilities_bgy",
  "capabilities_mun",
  "trust_bgy",
  "trust_mun")
state_imp <- imputeFun(all_nomat, state_vars)
prcomp_state <- prcomp(state_imp, scale = TRUE)

all_nomat$state_attitudes_index <- as.matrix(state_imp) %*% as.matrix(prcomp_state$rotation[, 1])*-1

# Estimate priming effects on index and individual attitudes
mod.prime.state_attitudes_index.id <- lm(state_attitudes_index ~ ipra_any, data = all_nomat)
mod.prime.intentions_bgy.id <- lm(intentions_bgy ~ ipra_any, data = all_nomat)
mod.prime.intentions_mun.id <- lm(intentions_mun ~ ipra_any, data = all_nomat)
mod.prime.capabilities_bgy.id <- lm(capabilities_bgy ~ ipra_any, data = all_nomat)
mod.prime.capabilities_mun.id <- lm(capabilities_mun ~ ipra_any, data = all_nomat)
mod.prime.trust_bgy.id <- lm(trust_bgy ~ ipra_any, data = all_nomat)
mod.prime.trust_mun.id <- lm(trust_mun ~ ipra_any, data = all_nomat)

# Create plot
pdf("state_attitudes_effects_id.pdf")
par(mfrow = c(1, 1))
modlist <- list(mod.prime.state_attitudes_index.id,
                mod.prime.intentions_bgy.id,
                mod.prime.intentions_mun.id,
                mod.prime.capabilities_bgy.id,
                mod.prime.capabilities_mun.id,
                mod.prime.trust_bgy.id,
                mod.prime.trust_mun.id)
res <- data.frame(matrix(NA, nrow = length(modlist), ncol = 3))
colnames(res) <- c("coef_id", "ci_lower_id", "ci_upper_id")
for(i in 1:length(modlist)){
  res$coef_id[i] <- summary(modlist[[i]])$coef[2, 1]
  se <- summary(modlist[[i]])$coef[2, 2]
  res$ci_lower_id[i] <- res$coef_id[i] - 1.96*se
  res$ci_upper_id[i] <- res$coef_id[i] + 1.96*se
}
modlabels <- c("index",
               "intentions \n bgy",
               "intentions \n mun",
               "capable \n bgy",
               "capable \n mun",
               "trust \n bgy",
               "trust \n mun")
plot(x = c(1:length(modlist)),
     y = res$coef_id,
     xlim = c(0.5, length(modlist) + 0.5),
     ylim = c(-0.7, 0.7),
     pch = 20,
     ylab = "coef",
     xlab = "",
     xaxt = 'n',
     main = "Experimental Effects on State Attitudes",
     cex.main = 0.8)
axis(side = 1, labels = modlabels,
     at = c(1:length(modlist)),
     cex.axis = 0.6)
segments(x0 = c(1:length(modlist)),
         y0 = res$ci_lower_id,
         y1 = res$ci_upper_id)
abline(h = 0, lty = 3)
dev.off()

# Figure A.7: Targey survey sample

# Province map
pdf("sample_map_prov_gray.pdf")
plot(prov_proj)
plot(prov_proj[which(prov_proj$LAYERBB %in% c("MINDORO_ORIENTAL", "MINDORO_OCCIDENTAL", "PALAWAN")), ], add = T, col = "grey36")
dev.off()

# Barangay map
pdf("sample_map_zoom_gray.pdf")
plot(prov_proj[which(prov_proj$LAYERBB %in% c("MINDORO_ORIENTAL", "MINDORO_OCCIDENTAL", "PALAWAN")), ], border = "white")
mimaropa_bgys <- bgy_map[which(bgy_map$prov_name_2000 %in% c("ORIENTAL MINDORO", "OCCIDENTAL MINDORO", "PALAWAN")), ]
plot(mimaropa_bgys, xlim = c(bbox(mimaropa_bgys)[1, 1] + 500000, bbox(mimaropa_bgys)[1, 2]), add = T)
plot(mimaropa_bgys[which(mimaropa_bgys$included_survey == 1), ], add = T, col = "grey36")
# add province labels
text("ORIENTAL MINDORO",
     x = coordinates(prov_proj)[which(prov_proj$LAYERBB == "MINDORO_ORIENTAL"), 1] - 180000,
     y = coordinates(prov_proj)[which(prov_proj$LAYERBB == "MINDORO_ORIENTAL"), 2] - 30000)
text("OCCIDENTAL MINDORO",
     x = coordinates(prov_proj)[which(prov_proj$LAYERBB == "MINDORO_OCCIDENTAL"), 1] + 70000,
     y = coordinates(prov_proj)[which(prov_proj$LAYERBB == "MINDORO_OCCIDENTAL"), 2] + 120000)
text("PALAWAN",
     x = coordinates(prov_proj)[which(prov_proj$LAYERBB == "PALAWAN"), 1] - 87000,
     y = coordinates(prov_proj)[which(prov_proj$LAYERBB == "PALAWAN"), 2])
dev.off()

